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PREFACE 


The outdoor propagation of sound remains an important topic of research. Some of the earliest recorded 
experiments in acoustics dealt with the propagation of sound. The reason for the continuing interest in sound 
propagation is that sound propagation is an aspect of many acoustic problems. In recent history, during 
the decade of the 70’s, outdoor sound propagation research was largely driven by aircraft noise certification 
issues. Propagation distances of interest were typically on the order of a mile. The effects of finite impedance 
boundaries, ground effects, were identified as important to the problem, and much theoretical and experimental 
work was done on ground effects. Today, propagation distances of interest are an order of magnitude larger. 
Propagation problems of interest include refraction due to speed of sound gradients and scattering due to 
turbulence. Applications of long-range sound propagation technology range from en route aircraft noise to 
the acoustic detection of aircraft. In 1981, the University of Mississippi and the Open University of England 
co-sponsored a symposium which dealt with issues of particular interest to outdoor, long-range propagation. 
Approximately every 2 years since the first, the University of Mississippi and the Open University of England 
have co-sponsored with a third institution a similar symposium on long-range sound propagation. The Fourth 
International Symposium on Long-Range Sound Propagation was held at the NASA Langley Research Center, 
Hampton, Va., on May 16-17, 1990. The purpose of the meeting was to exchange information on current 
research, identify areas needing additional work, and coordinate activities as much as possible. The list of 
attendees which follows includes representatives from most groups with active research programs in the area. 

The meeting was divided into three sessions: ground effects on propagation, infrasound propagation, and 
meteorological effects on sound propagation. The symposium ended with an open discussion and plans for a 
future meeting. This report consists of a list of attendees with addresses, a meeting agenda, and a compilation 
of the presentations made at the symposium. 

The hosts would like to express their appreciation to the participants for attending and for sharing their 
knowledge and expertise. 
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LONG-RANGE SOUND PROPAGATION - 
A REVIEW OF SOME EXPERIMENTAL DATA 

Louis C. Sutherland 
Consultant in Acoustics 
27803 Longhill Dr. 

Rancho Palos Verdes, CA 90274 


SUMMARY 

Three experimental studies of long range sound propagation carried out or sponsored in the past by 
NASA are briefly reviewed to provide a partial prospective for some of the analytical studies presented in 
this symposium. The three studies reviewed cover (1) a unique test of two large rocket engines conducted 
in such a way as to provide an indication of possible atmospheric scattering loss from a large low- 
frequency directive sound source, (2) a year-long measurement of low frequency sound propagation 
which clearly demonstrated the dominant influence of the vertical gradient in the vector sound velocity 
towards the receiver in defining excess sound attenuation due to refraction, and (3), a series of excess 
ground attenuation measurements over grass and asphalt surfaces replicated several times under very 
similar inversion weather conditions. 


INTRODUCTION 

Experimental data on long range sound propagation sound from three unique programs carried out 
over the last 25 years that were conducted or sponsored by NASA can provide a useful background for 
some of the analytical models treated in this symposium. These measurement programs are very briefly 
reviewed here to insure that the existence of these data may be more widely known to researchers in the 
field of long range sound propagation. The sources of the data are identified for the reader who may wish 
to pursue the information in more detail. 

EXPERIMENTAL DATA ON PROPAGATION OF 
LOW FREQUENCY ROCKET NOISE AT LONG RANGES. 

On March 24, 1964 at approximately 1340 CST, the NASA George C. Marshall Space Flight 
Center, in Huntsville, Alabama conducted a static test firing of a Saturn S-I first stage rocket booster on a 
test stand for which the deflected exhaust blast was directed due north. This rocket consists of a cluster of 
eight engines with a total thrust of about 1.5 million lbs. Seven minutes later, a static test of a Saturn F-I 
rocket engine (a single chamber rocket engine with the same total thrust), was conducted on the same basic 
test stand but with the deflected exhaust blast directed due south. Major results of acoustic measurements 
conducted out to a distance of 15 Km along a line of microphone stations on a 45° azimuth line from the 
test stand towards the city of Huntsville, as shown in Figure 1, were reported by Tedrick. 1 However, 
most of the detailed results presented here are contained in an internal NASA Memo. 2 Also shown in 
Figure 1 are the vertical sound velocity profiles measured in this direction at the time of each test firing and 
the resulting calculated sound ray paths in this same direction. The sound velocity profiles differ slightly 
in the first 2 Km but the resulting ray paths differ significantly. Based on a comparison of the ray paths 
for the two firings, one would expect to see a greater refraction loss for the second test due to the greater 
upward refraction of the sound ray for this test. As will be shown, precisely the opposite condition 
prevailed. 

Not shown here are the same type of sound profiles and ray paths for a 226° azimuth direction - 
essentially 180° from those shown in Figure 1. The results were very similar - minor differences in sound 
profiles and a ray path for the second test showing more upward refraction in this direction than for the 
first test — again suggesting a greater refraction loss for the second test. 



Although the two rocket boosters have a very different geometry, the resultant total sound power 
levels and spectra are very similar 1 and, as shown in Figure 2, the directivities for the overall sound 
pressure level at a distance of 1000 ft from the engines are very similar when the different direction of the 
exhaust blast for the two tests is recognized. In the direction of the microphone positions, the overall 
sound levels of the two rocket engines differ by about 12 dB at a 1000 ft radius. Figure 3 shows the 
values of excess attenuation in octave bands, including any air absorption, for the S-I test, as a function of 
octave band center frequency with distance as a parameter. It was convenient, for this plot, to use 1 .6 Km 
as a reference distance for evaluating excess attenuation. The data show, roughly, the expected trend of 
increasing excess attenuation with distance and frequency. Figure 4 shows the same data for the S-I test 
re-plotted as a function of distance where the values of excess attenuation have been averaged over pairs of 
adjacent octave bands to simplify the data presentation. Figure 5 shows the same information for the F-I 
test. 

However, it is not the purpose of this review to examine the absolute values for the excess attenuation for 
each test but rather examine the difference in excess attenuation between the two tests. This is shown in 
Figure 6 in terms of the excess attenuation for the S-I test (i.e., maximum lobe of noise along the 
measurement direction towards Huntsville) minus the excess attenuation along the same line, for the F-I 
test (i.e., maximum lobe of noise in opposite direction). 

The excess attenuation along this same path decreased between the two tests, conducted only 7 minutes 
apart. This decrease is most significant for a distance of 9 Km and is more dependent upon frequency at 
this distance than at any other point. This decrease in excess attenuation could be attributed to a change in 
sound refraction between the two tests. However, as suggested by the sound velocity profiles and 
calculated ray paths in Figure 1, this effect would have been expected to be just the opposite from what 
was observed - i.e., an increase in excess attenuation due to the expected increase in refraction loss for the 
second test. An alternative hypothesis is that the decrease in excess attenuation could be attributed to the 
effect of scattering by atmospheric turbulence. This scattering would tend to increase the apparent excess 
attenuation in the measurement direction for the first test (i.e., remove energy from the main sound lobe in 
this direction) and decrease the excess attenuation for the second test by adding back-scattered energy to 
the weaker lobe in this direction. 

This hypothesis, admittedly not proven, is consistent with the observations and with theoretical 
predictions. Further research is needed to more fully evaluate and experimentally validate sound 
attenuation by atmospheric turbulence. Practical applications include definition of correct excess 
attenuation models for the directive sound fields of jet aircraft and long range warning sirens. 

LONG-TERM MEASUREMENT OF EXCESS ATTENUATION 
WITH REFRACTION 

The second sound attenuation program was conducted at the NASA Mississippi Test Range over a 
one year period by Tedrick and Polly. ^ The program utilized the pure tone siren/hom sound source 
system shown in Figure 7 mounted on a 60 ft. tower to propagate pure tone signals at 40, 80, 120 and 160 
Hz at distances up to 3 Km over a flat terrain heavily covered with a deciduous rain forest. Over 29,000 
excess attenuation measurements were made over the one year test period. The results were correlated 
with the vertical gradient of vector sound velocity from the source to the receiver as measured over the first 
300 meters above the ground. Typical results for two distances are shown in Figure 8 in terms of the 
excess attenuation at 160 Hz as a function of this sound velocity gradient. As for all of the frequencies and 
distances measured, the data collapsed in the form illustrated. At any given frequency and distance, the 
mean excess attenuation was essentially constant when the sound velocity gradient was equal to, or greater 
than zero and decreased approximately linearly as the gradient decreased below zero. 

The mean excess attenuation, A 0 for sound velocity gradients equal or greater than zero varied 
linearly with distance and systematically with frequency as shown on Figure 9 which is taken from Ref. 5. 
Although the excess attenuation includes air absorption, the latter is a relatively small part of the observed 
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excess attenuation which is believed to be predominantly ground attenuation. Note that the intercept value 
of Aq for zero distance is roughly proportional to frequency but the rate of increase with distance increases 
only slightly with frequency. 

For negative sound velocity gradients, Tedrick and Polly showed that the slope of the plot of 
excess attenuation versus sound velocity gradient increased linearly with distance and approximately 
linearly with frequency (see Figure 10). 

While the above presents a very simplified definition of the data trends, it has substantial face 
validity on the basis of the very large number of measurements involved and should provide useful 
benchmarks for comparison with the latest theoretical models for ground attenuation in the presence of 
refraction. 

Another result from this long term test program was the determination of the statistical distribution 
in the magnitude of focusing amplification (i.e., excess attenuation which is positive) corresponding to 
sound attenuation less than inverse square spreading loss. While very likely a site-specific statistic, the 
distribution data shown in Figure 11, developed from tabular data in Ref. 5, shows that this focusing 
anomaly increases with distance for values of the anomaly less than about 15 dB. Note that in this case, 
the data cover a much longer propagation range and indicate that, on rare occasions, anomalous increases 
in level above that predicted by spherical spreading loss of up to 30 dB were observed. 

GROUND ATTENUATION MEASUREMENTS FOR INVERSION CONDITIONS 
OVER GRASS AND ASPHALT SURFACES. 

The final test program mentioned here was sponsored by NASA and is fully described in Ref. 6. 
Copies of the full report may be available through NASA, Langley. The program involved the 
measurement of ground attenuation over asphalt and grass surfaces on, or next to, an aircraft runway at 
NASA's Wallops Island facility. The tests were conducted with an elevated loudspeaker source located at 
2.5, 5, and 10 meters above each of the surfaces. For most of the tests, the weather conditions 
corresponded to a mild inversion condition that was replicated several times for each measurement source 
elevation/ground surface condition. The basic test geometry and microphone array employed is illustrated 
in Figure 12. Note that at one distance (225 meters), microphones were located essentially at the ground 
surface, and at 1.2 and 10 meters. At 450 meters, microphones were located at 1.2 and 10 meters. (Note 
that for the tests over grass, a small strip of asphalt existed along the "grass" path between the 450 and 675 
m positions.) 

Along with the excess attenuation measurements, the mean weather conditions were evaluated 
extensively with meteorological instrumentation on 7 and 10 meter towers and a captive weather balloon 
repeatedly raised to and lowered from a height of 100 m. For the sake of brevity, only a small fraction of 
the available excess attenuation data are shown here in Figure 13. The figure shows, for two distances, 
the two surfaces and three source heights, the arithmetically averaged excess attenuation for one third 
octave bands of noise from 50 to at least 3200 Hz for the four to six replications of nominally very similar 
inversion conditions. Each excess attenuation measurement was based on an energy average of sound 
levels oyer a 1 5 second period. The standard deviation of the excess attenuation values over the four to six 
replications for each measurement condition and frequency was normally much less than 1.5 dB. 

The results show the characteristic increase in excess attenuation due to ground absorption at frequencies 
in the range of 125 to 630 Hz depending on the surface and measurement distance. The excess attenuation 
data are augmented by some very limited measurements of surface impedance employing the simple 
technique developed by Piercy and Embleton. 7 Thus, these data provide another, and, in some aspects 
more complete, set of measurements of ground attenuation in the presence of documented refraction 
conditions than had been available previously. They offer a useful set of measurements for comparison 
with corresponding theoretical models. 
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CONCLUSIONS 


Results from three different NASA conducted or NASA sponsored tests of long range sound 
propagation have been very briefly reviewed. The objective has been to identify these unique sources of 
data, two of which are over 25 years old, for the benefit of modelers of long range sound propagation who 
may not be aware of their existence. They offer potentially useful data sets for comparison with theoretical 
models for the evaluation, respectively of: scattering attenuation by atmospheric turbulence, long range 
ground propagation under a wide range of defined refraction conditions, and ground attenuation over two 
surfaces for nearly identical mild inversion conditions. As further advances are made in theoretical 
models, new and more sophisticated measurements will be required to validate the theory. 
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Figure 9. Excess Attenuation for Positive Sound Velocity Gradient 



Figure 10. Slope of Excess Attenuation for Negative Sound Velocity Gradient 
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Figure 12. Ground Attenuation Measurements at Wallops Island (Taken from Sutherland 
and Brown, NASA CR-3435, 1981) 
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ABSTRACT 


The propagation of grazing incidence plane waves along a finite impedance boundary is 
investigated. A solution of the semi-infinite problem, where a harmonic motion, parallel to the 
boundary, is imposed along a line perpendicular to the boundary, is obtained. This solution consists 
of quasiplane waves, waves moving parallel to the boundary with amplitude and phase variations 
perpendicular to the boundary. Several approximations to the full solution are considered. 


INTRODUCTION 


Mathematical modeling of the propagation and reflection of harmonic plane waves above a 
finite impedance plane surface is a fundamental topic in acoustics. In the case where the angle 
between the normal to the wavefront and the surface is not zero an analytic solution is very easy to 
obtain. This solution consists of the incident plane wave propagating toward the surface plus a 
reflected plane wave propagating away from the surface at the same magnitude of the angle between its 
normal and the surface as the incident wave. The amplitude of the reflected wave is given by a 
reflection coefficient that is expressed in terms of incident angle and the specific impedance of the 
surface. However at zero incident angle (the wave normal parallel to the surface), complete cancellation 
of the incident and reflect waves occurs in this model and a zero solution results. Most acoustic texts 
claim that this situation is not possible [1-3]. 

McAninch[4] recently has investigated a related situation where a plane wave source is 
generating waves that would move parallel to a surface if its impedance was infinite but where the 
surface impedance isnot infinite quasiplane waves result McAninch's investigation, however, uses the 
parabolic approximation where only waves traveling in one direction are allowed. This paper 
approaches the same problem without the assumption of parabolic approximation. 


FORMULATION OF THE PROBLEM 


The governing acoustic wave equation for harmonic waves can be put in the form 

( V 2 + k 2 ) <J>=0 (1) 

where the time dependent part of the potential, e- i(0t , has been separated from the spatial part of the 
potential, 4>(x,y). When an impedance boundary exists, the solution of equation (1) must also satisfy 
the boundary condition 

4> y + Y<1> = 0 (2) 

on y=0. Here the subscript y indicates a partial derivative with respect to y. 
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For a unique solution, some extra constraints must be introduced. One is to assume that <|) will 
not be affected by the ground impedance as y approaches infinity, the second is that the acoustic 
pressure at x=0 is given by 


<K0,y) = 1 


(3) 


EXACT SOLUTION 


We assume that the solution of (1) has the form of 

([) = e ikx + f(x,y) (4) 

where e ikx can be considered as a solution without the boundary condition given by (2). Substituting 
(4) into (1),(2) and (3) we get a new governing equation and set of boundary conditions 

( V 2 + k 2 ) f = 0 (5) 

f y (x,0) + yf (x,0) = - ye 1 k x (6) 

f (0,y) = 0 (7) 

and 

lim f (x,y) = 0 (8) 

y— 

Equation (8) results from the first uniqueness condition listed above. 

The sine transform, 

oo 

F (X,y) = J f (x,y) sin (kx) dx (9) 

o 

is equivalent to the Fourier transform of an even function and will be applied here. The inverse 
transform is given by 


f (x,y) = ~ J F (x,y) sin (kx) dk 
0 


Applying (9) to (5) and (6) we have 


and 


3 2 F 

3y 2 


( k 2 - k 2 ) F = 0 


F (k, 0) + yF (k,0) = - y— r^— r 

7 1 -k 2 


( 10 ) 


(ID 

( 12 ) 
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Here, it is assumed that k is a complex number with a very small positive imaginary part. 
Solving (11) and making use of the given boundary conditions, yields 


F (X,y) = A (X,0) e my 

(13) 

A(K 0 )~- 

( X. - k ) ( y- m ) 

(14) 

where m = V(X 2 -k 2 ) . Since the solution is required to remain finite, ReV(X 2 -k 2 ) >0. Substituting (13) 
and (14) into (10) yields the inverse transform of F(X,y) as 

oo 

f (x,y) = - -rp f — — — e my sin (X,x) dX. 

" i(X 2 - k 2 )(y-m) 

(15) 

For convenience, substitute the identity 


. J Xx i Xx 

sin (Xx) = — — ^ 

2 1 

into (15), yielding 

(16) 

f( x,y ) = - 1 X (I ,- I 2 ) 

(17) 

where 


oo 

r 1 ’ my 

•l = J“T! eiX%dX 

o (*> ' k )(7-m) 

(18) 

oo 

r 'l _ m y 

>2 = [—5-“ 

o (^ - k )(y-m) 

(19) 

In order to evaluate the above two integrals, introduce the complex variable A = 

X+is and define the 


contour integrals 



Ae 


-My 


( A 2 - k 2 ) ( y- M) 


Ae 


-My 


( A 2 - k 2 ) ( y- M) 


e iAx dA 
e 1 Ax dA 


( 20 ) 


( 21 ) 


First evaluate the integral Iq where the contour is shown in Figure 1 along with the branch lines 
which extend from the imaginary axis to the points A = ± k. The value of this contour integral is 
determined by the residue within the contour. Writing 


Y= a + i |3 


( 22 ) 
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it is clear a pole exists within the contour only when a > 0 (since Re^A i 2 - k 2 > 0, y - VA 2 -k 2 =0 only 
when Re(y) = a > 0), and this pole is at A = \(k 2 + y 2 ). It is easy to determine the residue at this pole 


to be 


Res 


(Vk 2 +f )=-£ 


yy +i 7 k 2 +-y ? x 


and 


*1 " *C1 ' ^C2 ' ^C3 ' *C4 ' ^C5 ' *C6 + 


2 7t i Res ( -J k 2 + y 2 ) a > 0 
0 a < 0 


when R — > co, l C2 will vanish, while 

OO 

T C3 = ' J 


-i 7 s 2 + k 2 


y e sx ds 


o ( s 2 + k 2 ) ( y- i 7 s 2 + k 2 


(23) 


(24) 


(25) 


C4 


K. 

"J 


-i 7k 2 - X 2 y j Xx 


2 ' 2 ^-- iTi 2 ^ 2 


e' dX 


o ( X - k 6 ) ( y- i V k' - X ) 


t - . ILL J k x 
X5 y e 


I 


C6 


K 

-J' 


i 


‘ 2 ' 2 ivV 


e 


o (X. -k ) (y+iVk -X. ) 

Substituting the above integrals into (24), yields 


K 




.2 \ 2 
k - A, y 


-i 7k 


k 2 -\ 2 y 

6 ^ e‘ Xx dX 


J(?i 2 -k 2 ) [(y+i7k 2 -^ 2 ) ( y- i Vk 2 - X. 2 ) 


i 7 s 2 + k 2 


oo 

I — 

o ( s 2 + k 2 ) ( y- i v s 2 + k 2 


e' s x ds + 2-ZL e ikx 


+ 


i 2 n - yy J -/kV^ 

Y 

0 


a > 0 
a < 0 


(26) 


(27) 


(28) 


(29) 
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The value of integral I 2 is much easier to evaluate. We choose a contour in the fourth quadrant, since 
there is no pole within the contour. I 2 can then be written as 



s e 


ij 


s 2 +k 2 y 


( s 2 +k 2 ) ( Y+iTs 2 


+ ) 


e' sx ds 


(30) 


Here it must be recalled that there is a branch line along the imaginary axis. Subtracting I, from I, 
yields 



- i V s 2 + 


k 2 y 


i 777 


k 2 y 


00 

f s 

o ( s 2 + k 2 ) [ ( y- i V s 2 + k 2 ) ( y+ i 7 s 2 + k 2 ) 


>e sx ds 


+ 




i 271 -yy i •/ k 2 +y 2 x 

y c c 
0 


a > 0 
a < 0 


(31) 


By substitution of A.=Vk 2 -t 2 and s=-W k 2 -t 2 , the above two integrals can be combined into one. The 
final result is 


where 


and 


f(x,y) 


- e‘ k x + P - K 
-e ikx -K 


a > 0 
a < 0 


P = 2 e 


- yy 



X 


(32) 


(33) 


K= 2jr f ( t Cos (ty) - YSin (ty) ) ^ V k 2 . t z x dt 
71 g t(f+t 2 ) 


(34) 


P is called the surface wave, and it both decays with increasing height y, and also decays with the 

distance x due to the imaginary part of Y 
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ASYMPTOTIC EXPANSION VALID FOR SMALL RECEIVER HEIGHTS 


1. Soft boundary case 

Integral K can be asymptotically evaluated for large x using the saddle-point method. This 
method is discussed by Morse and Feshbach[5] and will not be discussed here. Actually we can use 
some of conclusions from Wenzel [6] since we have the same factor Vk 2 -t 2 as occurred there. 

The steepest-descent path has been shown in [6] to be given by 

T = t + i s (35) 

where 

(36) 


and t > 0. Again using the residue theorem, integral K can be transformed into the integral L. Note 
that i y is not in the region of concern since |3>0. Thus 



K = <! 


f P + L 


- i ye D 

- i Y £ D 


(37) 


where D is the region between positive real axis and curve s=-t (1+t 2 / k 2 ) 1/2 , P is the surface wave 
given in (33) and 


2 Y f T Cos (Ty) - y Sin (Ty) x dT 

71 sit (f+T 2 )T 


(38) 


Substituting (37) into (32), we have 



P + L 
L 


- i y e T 

- i Y <£ r 


(39) 


The region T in y plane is bounded by the curve p > 0, 0 < a=p(l+ (3 2 / k 2 ) 1/2 . The region T is 
called the surface wave region in the far field (shown in Figure 2), which is same as that of reference 
[6]. When ye T, we can easily show ReV k 2 + y 2 > k,that means, if the surface wave exists,its 
propagation speed is less than the speed of sound in free space. It is also found that Im "v k 2 + y 2 
has a close relationship to the quantity (a P / k) so that a large imaginary part of y and low frequency 
of the source can make the surface wave decay very quickly. 

Expanding each expression in (38) around the saddle point T=0 and integrating each term, we 
get 
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( 40 ) 


=i /jO 

Y V in 


i lc T 3 y 3 9 

Jkx<!(l-Yy) + -LL (1- Yy+ 1JL.JL^) + 0(X- 2 ) 
I x 


This asymptotic expansion is not uniformly valid. The conditions for its validity are 


and 



(41) 

(42) 


Extremely small l Y I will not satisfy conditions (41) and (42), so another asymptotic method has to be 
developed. 

The total solution under the condition of large x can be obtained by substituting (39) into (4) 
yielding 


f P-L 



y € r 

Y e r 


(43) 


where P and L are given in (33) and (40). If we neglect the surface wave, we can get an explicit 
equation for the wave above ground in the far field as 


(j> = 


lYl 


7 ? 


((l-ay) 2 + (Py) 2 ) 


i(kx+^- +0) 
_ 4 


(44) 


where 0 = - arctan (p y / (1 - a y)). Furthermore, if the receiver is on the ground, the above 
expression can be written as 


20 Log <J>= 20 Log a - 10 Log x (45) 

with a = (l/lY0(2k/7i) 1/2 . This result shows that the acoustic pressure level drops lOdB when the 
distance increases 10 times or 3 dB per doubling of distance. 

2. Hard boundary case 

As mentioned before, the asymptotic expansion given in (40) is not uniformly valid in y, with 
the method failing for small lyl. An alternative method is developed in this section which is valid in 
the small lyl case. The method is almost the same as that used in evaluating L except the factor 1 / 
(Y 2 + T 2 ) in (38) will not be expanded. After changing variables (38) becomes 


L = 


Y ( 1 - YY ) 

7t i 


Making use of the formula 



A kx 


I< 


X “jp 

2 i k 



(46) 
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( 47 ) 


yry e~ z2 erfc ( -iz ) Im (z) > 0 

-~-t z ( erfc ( -iz ) - 2 ) lm(z)<0 

2 1 z 

and neglecting the terms of order y 2 , yields 

-e ikx [l+y(-l r /IT _ y)] Im(-?-)>0 

Jti v 2 k yi 

(48) 

. e i kx [i + y(^_ /IT -y)-2(l-yy)] Im(-|-)<0 

V TC V 2 K yj 


The conditions Im (y/ Vi ) > 0 and Im (y/ Vi ) < 0 can be identified as a < (3 and a > (3 respectively, 
a = (3 is the line which divides these two regions in y plane. This is exactly the bounding curve a = [3 
(1 + P 2 /k 2 ) _1/2 obtained previously provided that lyl — » 0. Recognizing this relation, we substitute (48) 
into (43) and rewriting surface wave approximately as P = 2 (1 - y y) e ik \ finally get the total field 
expression as 





the condition for the validity of the above expansion is 

lY|y«l 

and 

lir | V ! « 1 

although x can't be small because of the nature of the saddle point method. 


(49) 


(50) 


(51) 


ANOTHER ASYMPTOTIC EXPANSION VALID FOR LARGE RECEIVER HEIGHT 


The asymptotic expansions obtained above have their limitations in application. For example, 
they require the receiver's location to be near the ground. In this section we will derive a asymptotic 
expansion which is valid for large R = Vx 2 +y 2 (except for small y) . The idea is similar to that of 
Chien and Soroka [7]. 


By using the identity 

f(x,y) 


sin (Xx)=(e‘^ x -e' i ^ x )/2i and transformation 
_ y f Tan ( z ) e ik ^ y Cos(z)+x 

in J y + i k Cos ( z ) 


X=k sin (z), (15) becomes 
s,n<l)) dz (52) 
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The contour C is shown in Figure 3. In order to get an expansion in terms of the variable R = Vx 2 + y 2 , 
we transform the Cartesian coordinate system into the polar coordinate system by 

x = R Sin 0 (53) 

and 

y = R Cos 0 (54) 

Substituting into (52), yields 


f (x,y) = -X- f — ^ 

1 71 J Y+ 1 


Tan (z) 


y + i Cos ( z ) 


ikRCos(z-G) j 

e dz 


(55) 


the saddle point for the function i k R Cos(z-G) is at z = 0 and the path of steepest descent is found to 
be given by 


Cos ( u - 0 ) Cosh ( v ) = 1 (56) 

where z = u + i v, by considering Im (i k R Cos (z - 0))=lm (i k R). This path, denoted as C' is 
shown in Figure (3). Deforming contour C into C', adding the possible poles (Cos z = i y / k ), we 
have 


f (x,y) = Q + H ( - Re ( 1 - ^ Cbs 0 - 

Jv 



Sin 0 ) ) P 


(57) 


where Q is defined by (54) but with the contour C changed to C', H is Heavyside step function and P 
is the surface wave given in (33). The condition for the existence of pole is explained in reference [7], 
and will not be repeated here. In the limit of 0 approaching 7t/2 the condition for the existence of the 
pole in the present case is equivalent to the condition for the existence of the pole in (39). 

Q can be evaluated asymptotically with a method similar to that used in evaluating L, i.e. to 
expand each term around the saddle point 0 and then integrate them with suitable transformation of the 
variable . The result is 


q= /24— — — 

V 1 7t k R i ( y+ i 


yTan 0 ikR 
c 


1 +T 


1 


ikR 


1 + 


Cos 0 ) 
i k Cos 0 


i k 


2 ( y+ i k Cos 0 ) ( y + i k Cos 0 ) Cos 0 


The conditions for the validity of the above expansion are 

kR » 1 


and 


Cos 2 0 , ^ 


y + i Cos 0 R Cos 0 


« 1 


(58) 


(59) 

(60) 


21 



It is clear that 0 cannot be too close to Till. This limitation is complimentary to the asymptotic 
expansions obtained previously (for small y). Substituting into (4) yields 


iJ)=e ikx +Q + H(-Re(l - Cos 0 - 1 + “T Sin 0 )) P (61) 

In the limit R->oo Q and P will vanish , with the result that only the plane wave term remains. 

Equation (38) can be evaluated accurately by numerical methods as well as by asymptotic 
expansions. Calculations show that the results match quite well when y is small. Figure 4 gives the 
amplitude of acoustic pressure on the ground versus the distance to the receiver obtained by numerical 
integration and from (44). Figures 5a, b and c show the amplitude of acoustic pressure versus the 
receiver height for several receiver locations as obtained from the asymptotic expansions, (40) and 
(61). These figures are similar to the results obtained by McAninch [4], 


CONCLUSIONS 


The acoustic field of a plane wave at grazing incident to a finite impedance has been 
theoretically investigated. Exact numerical and asymptotic expansions are developed, which are very 
similar to those found by Wenzel [6] for a point source and by McAninch f4j using the parabolic 
approximation to the wave equation. When y is small, the incident wave is indeed canceled, but the 
result is not zero due to the existence of a surface wave and the wave denoted as L. Near the ground, 
the acoustic pressure decays as x' 1/2 (assuming the surface wave is neglected). The asymptotic 
expansion for large distance R shows that the acoustic pressure decays as R 1/2 when R— and 
when the receiver is not close to the surface only incident wave exists. 


REFERENCES 


[1] Dowling, A. P., and J. E. Ffowes Williams, Sound and Sources of Sound, (Ellis Horwood Ltd., 
Chichester, 1983), p. 83. 

[2] Morse, P. M., Vibration and Sound, (American Institute of Physics, Woodbury, N.Y., 1983), 
pp. 366-368. 

[3] Rudnick, I., "The Propagation of an Acoustic Wave along a Boundary." J. Acoust. Soc. Am., 
19 , 348-356, (1947). 

[4] McAninch, G. L., "Propagation of Quasiplane Waves Along an Impedance Boundary." AIAA 
26th Aerospace Sciences Meeting Paper 88-0179 (1974). 

[5] Morse, P. M and H. Feshbach, Methods of Theoretical Physics, (McGraw-Hill, New York, 1953), 
p. 441. 

[6] Wenzel, A. R., "Propagation of Waves along an Impedance Boundary." J. Acoust. Soc. Am., 
55,956-963 (1974). 

[7] Chien,C.F. and W.W. Soroka, "Sound Propagation Along an Impedance Plane." J. Sound and 
Vibration, 43 ,(1), 9-20, (1975). 




Figure 3. The u-v plane, showing the integration contour C and the steepest descent path C'. 
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Figure 4. Comparison between a numerical result and the asymptotic result, (44), for the parameters, 
k = 18.313 and y = 6.935 + i 19.015. ’ ‘ 
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ABSTRACT 

In atmospheric acoustics, the subject of surface waves has been an area of discussion for 
many years. The existence of an acoustic surface wave is now well established theoretically. 
The mathematical solution for spherical wave propagation above an impedance boundary 
includes the possibility of a contribution that possesses all the standard properties for a 
surface wave. Surface waves exist when the surface is sufficiently porous, relative to its 
acoustical resistance, that it can influence the airborne particle velocity near the surface and 
reduce the phase velocity of sound waves in air at the surface. This traps some of the 
sound energy in the air to remain near the surface as it propagates. Above porous grounds, 
the existence of surface waves has eluded direct experimental confirmation (pulse experiments 
have failed to show a separate arrival expected from the reduced phase speed) and indirect 
evidence for its existence has appeared contradictory. In PART I of this paper the 
experimental evidence for the existence of an acoustical surface wave above porous 
boundaries is reviewed. Recent measurements including pulse experiments will also be 

described. 

A few years ago the acoustic impedance of a grass-covered surface was measured in the 
frequency range 30 to 300 Hz. In PART II of this paper further measurements on the same 
site are discussed. These measurements include core samples, a shallow refractive survey to 
determine the seismic velocities, and measurements of the acoustic-to-seismic coupling 

coefficient. 


PART I 


INTRODUCTION 

In atmospheric acoustics, the subject of surface waves above porous grounds has been 
an area of discussion for many years. The existence of an acoustic surface wave is now well 

established theoretically. The mathematical solution for spherical wave propagation above an 

impedance boundary includes the possibility of a contribution that possesses all the standard 

properties for a surface wave. These include cylindrical spreading in the horizontal 

direction, exponential decay in amplitude with height above the ground, and a reduced phase 
speed. 


However, above natural porous ground surfaces, the existence of an acoustic surface 
wave has eluded direct experimental confirmation. Pulse experiments have failed to show a 
separate arrival from the direct pulse as expected from the reduced phase speed. Further, 
indirect evidence for its existence has appeared contradictory. 


PRECEDING PAGE BUNK NOT FILMED 
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The experimental evidence for surface waves has been mostly restricted to careful 
indoor measurements, using sources of continuous sound and model surfaces composed of a 
thin layer of porous material or comblike structures. The reduced phase speed and 
cylindrical spreading of the surface wave are expected to produce a total sound pressure level 
in excess of that which would be measured over an acoustically hard boundary. 

In this paper the experimental evidence for the existence of an acoustical surface wave 
above porous boundaries is reviewed. In addition, some recent measurements including pulse 
experiments will also be discussed. 


FIGURE 1 

At this point it is useful to distinguish between body waves and boundary waves. 
Acoustic waves propagating through the body of the fluid are referred to as body waves. 

The effect of boundaries upon these waves is secondary in that the existence of the waves is 

in no way tied to the presence of the boundaries. The role of boundaries is strictly 
extrinsic. On the other hand, boundary waves depend upon the existence of boundaries to 
support them and the role of the boundaries here is intrinsic. 

In atmospheric acoustics, the field from a point source above a porous ground is 
commonly described in terms of direct, reflected, ground, and surface waves. Obviously 

ground and surface waves are closely related but their fundamental origins differ, as does 
their behavior during propagation. Ground waves exist because curved wave fronts strike 

different parts of the ground at different angles of incidence and because the reflection 

coefficient of finite-impedance ground is also a function of angle of incidence. Ground 

waves exist unless the ground is infinitely hard or infinitely soft or unless the incident wave 
fronts are plane, that is, the source can be considered infinitely far away. Ground waves 
can exist in the absence of surface waves. 

Surface waves exist when the ground surface is sufficiently porous, relative to its 
acoustical resistance, that it can influence the airborne particle velocity near the surface and 
reduce the phase velocity of sound waves in air at the surface. In its simplest terms, the 

condition for its existence is when the imaginary component of the surface impedance is a 

spring-like reactance and is greater than the resistive component. This traps some of the 

sound energy in the air, regardless of the shape of the incident sound field, to remain near 
the surface as it propagates from the source to the receiver. Surface waves can exist in the 
absence of ground waves. The existence of a surface wave in the absence of wavefront 

curvature has been shown theoretically by McAninch and Myers (AIAA 1988). They 
demonstrate the presence of a surface wave in the solution for plane waves at grazing 

incidence to a finite impedance boundary. Further, Raspet and Baird (JASA 1989) have 
demonstrated that the surface wave can exist independent of the acoustic body wave in the 

half-space above the surface by examining the limit as the upper half-space becomes 

incompressible. 

The equation on the top part of this figure represents a particular representation for 
the total field above an impedance plane. The field is broken up into a direct wave, a 
perfect reflected wave, a diffracted wave that accounts for the phase change on reflection 
and the effects of the spherical wave fronts, and a surface wave. The surface wave exists if 
Im(Z) > Re(Z) and is zero otherwise. The surface wave is characterized by cylindrical 

spreading in the horizontal plane, exponential decay with increasing height above the ground, 
and a reduced phased speed v < c. 

Theory which predicts the acoustical characteristics of rigid porous materials in terms 
of their microstructure indicates that the resistive and reactive components of the surface 
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impedance are equal in the case of a homogeneous porous ground (Attenborough, JSV 1985). 
Therefore, no surface wave can exist above such grounds. On the other hand, if the 
microstructural properties of the ground vary with depth (such as a varying porosity), the 
reactive component of the impedance exceeds the resistive component and the surface wave 
can exist. 

A specific example of a surface whose reactive component of impedance exceeds the 
resistive component is a thin porous layer above an acoustically hard backing. We note that 
in the case of a ground where the porosity varies with depth at a rate a, the impedance is 
equivalent to the impedance of a porous layer with an effective thickness equal to 2/a 
(Donato, JASA 1977). 


FIGURE 2 

The consequence and origin of the reduced phase speed of the surface wave are 
illustrated in this figure. Far from the ground, there is horizontal particle motion associated 
with the propagating body wave, as shown in A. Due to the alternating compression and 
rarefaction cycles, the air molecules at the ground are entrained in vertical particle motion 
as shown in C. Just above the surface of the ground in the fluid, the resulting particle 
motion is therefore elliptical, as shown in B. 

The elliptical particle motion results in a reduced phase speed and the resulting lag 
causes the wavefronts to be "bent" towards the ground, giving rise to enhanced sound energy 
close to the surface. The increased sound energy associated with the surface wave close to 
the ground is at the expense of less sound energy at heights above the surface wave 
thickness. This will be illustrated in some of the following figures. 

FIGURE 3 

This figure shows experimental evidence measured outdoors over natural ground 
surfaces. The points in (a) are measurements obtained by Rasmussen above grass covered 
ground. The sound pressure levels in this figure, and all of the following figures, are plotted 
relative to free field. Hence, these results suggest sound pressure levels in excess of the +6 
dB expected at lower frequencies. The solid curve is the best prediction that can be 
achieved by assuming the ground to be a semi-infinite half plane. Rasmussen calculated the 
dashed curve by assuming a porous layer 0.01 m thick. Equivalently, the same result can be 
obtained by assuming a ground with its porosity varying with dept at a rate given by a = 
2/0.01 = 200 m~l. This is a more likely physical model for natural ground surfaces (Donato 
JASA 1977). 

We note that the behavior of Rasmussen’s measurements is consistent with the 
behavior of the classic measurements of Parkin and Scholes (JSV 1965) of the propagation of 
jet engine noise above grass covered airport ground. 

In (b), the points were measured above a well defined layer (8 cm) of snow above 
frozen ground. The dashed curve was calculated by assuming a layer of snow infinitely 
thick. The solid curve accounts for the layer. Although the measurements show the enhanced 
dip that is predicted around 300 Hz (Chien and Soroka, JSV 1975), the behavior of the 
measurements at the lower frequencies indicate that the surface wave is absent. These 
results have contributed to the controversy concerning the existence of surface waves above 
natural ground surfaces. It has been suggested (Attenborough, JSV 1988) that the situation is 

complicated by the existence of seismic quarter-wavelength resonances in the low frequency 
range as a result of the elasticity of the porous surface layer. 
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FIGURE 4 


The short dashed curve on this slide is the sound pressure levels predicted for 

propagation at grazing incidence above an infinitely thick surface of porous felt. The 

propagation distance is 2 m. There is no surface wave and this curve represents the ground 

wave. The open squares are measurements obtained above a thick layer of felt. 

The upper two curves are calculated from different versions of the same theory in the 
case of a layer of felt of thickness 0.003 m. In this case the surface wave $ s exists. The 

difference between the two curves is attributed to numerical precision and is not significant 

for the discussion here. 

The solid points are measurements made by Thomasson above a layer of felt. The 

open circles are our own measurements and confirm the results of Thomasson. Both theory 
and experiment clearly indicate sound pressure levels in excess of the +6 dB expected from 

inverse square law above a perfectly rigid ground. 

FIGURE 5 

The open points are measurements made as a function of height above the same layer 

of felt and shown for two frequencies. The solid points were obtained by Thomasson for the 
same two frequencies. The solid and broken curves are the predictions calculated from two 

versions of the same theory. The broken curves are the predictions in the case of an 

infinitely thick layer. 

The dotted lines drawn at +6 dB show the levels expected in the case of a perfectly 

rigid ground. Both theory and measurements show the existence of the enhanced sound 

levels at heights below 10 cm resulting from the existence of the surface wave. In addition, 

the slightly reduced levels above about 10 cm, especially at 2 kHz, is observed. 

FIGURE 6 

i 

In this figure, the porous layer is replaced by a comblike surface consisting of 

overhead lighting panels (Donato, JASA 1978). The panels are molded plastic: there is a 

square array of solid ribs at 1.13 cm spacing; the sheet is 2.26 cm thick, open on top and 
bottom surfaces. The sheet is laid on a hard floor. 

Results of measurements are shown for two frequencies and two distances of j 

propagation. The solid points clearly show significantly enhanced sound levels close to the 

surface, especially at 800 Hz, and the expected reduced level at higher heights. 

The open points are the results above a rigid surface and the solid lines are drawn at 

+6 dB. 


FIGURE 7 

These results are similar to the ones on the previous slide but the first four meters of 
the propagation path are acoustically rigid while the remainder consist of the comblike surface. 

The solid points to the left are measurements made above the rigid surface. The open 
points on the right were measured 5 m from the source, hence after 1 m of propagation 
above the ceiling panels. 

The behavior of the results at 5 m suggest that a surface wave has developed over the 
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1 m of panel. We note that the panels are located about 12 wavelengths from the source. 

Therefore the surface wavelike behavior is exhibited when the curvature of the wavelength is 
significantly reduced. This is consistent with the theory of McAninch and Myers. 

FIGURE 8 

The solid points on the top part of this figure are the results measured at grazing 
incidence above the comblike surface as a function of frequency for a distance of 1 m. The 
behavior of these results is identical to those measured above the layer of felt. 

The solid curve on the bottom part of the figure shows the predicted surface wave 
velocity, v (Brekhovskikh, Sov. Phys. Acoust. 1959). The straight line at about 340 m/s 

indicates the speed of the body wave in air. Beyond about 1.5 kHz there is a sufficient 

difference between the surface wave velocity v and the body wave velocity c, that it should 
be possible to observe the surface as a separate arrival using a short pulse of sound 
propagating over a distance of a few meters. 


FIGURE 9 

The traces shown here are of a 2.1 kHz tone burst measured after propagation above 
the comblike surface at various distances up to 1.5 m. The arrow immediately below the 
last three traces indicates the arrival of the surface wave relative to the body wave 
predicted from the solid curve on the previous slide. 

The observed behavior of the measured pulses as a function of distance is not 
inconsistent with expectations. In the absence of a surface wave all the traces would have 
the appearance of the top trace. 


FIGURE 10 

This figure shows the traces at different receiver heights for three distances of 
propagation (the source is on the ground). At a distance of 0.1 m, the surface wave has not 
yet had time to develop and the trace does not change with height. 

At the other two distances, the exponential decay of the second arrival as a function 
of height is clearly illustrated and is indicative of a surface wave. 


PART II 


INTRODUCTION 

A few years ago the acoustic impedance of a grass-covered surface was measured 
(Daigle and Stinson, JASA 1987) in the frequency range 30 to 300 Hz by measuring the 
pressure, phase and phase-gradient in the sound field along a vertical line directly below a 
loudspeaker suspended some 7 m above the surface. Recent core samples showed that this 
ground consisted of a layer of silt of uniform texture and almost constant thickness (1.6 +/- 
0.3 m) over bedrock — a ground structure of ideal simplicity for acoustical study. Seismic 
velocity measurements were consistent with this simple structure, and indicated a layer 
thickness (1.9 +/- 0.3 m) reasonably in agreement with the core sample. 

The calculated quarter- wavelength-layer- thickness frequency is then about 45 Hz. 
Direct measurement of the acoustic-to-seismic coupling coefficient at normal incidence shows 
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maxima in the admittance of the surface at about 50 and 135 Hz. (Several other maxima 

exist at apparently unrelated frequencies.) At oblique angles of incidence the admittance 
spectrum is of similar shape but shifts upwards in frequency by about 10%. 

A number of minima in the admittance spectrum are also present and should correspond 

with maxima in the acoustic reflection coefficient; however, the correspondence was found to 

be poor. Probable explanations of the discrepancies could be that the ground exhibits in 

reality a more complex structure than our current understanding allows or that different 

measurements were over slightly different areas of the ground and detected different 

thicknesses of the supposedly constant thickness silt layer. 

FIGURE 11 

This figure illustrates the original measurements. A pure tone is radiated spherically 
from a loudspeaker suspended resiliently from a support. Wavefronts are reflected at the 
ground surface and interfere with the incoming waves to produce an interference field. Two 

closely spaced microphones were moved together along a track that was perpendicular to the 
surface and directly below the source. By comparing the signals from the two microphones 
with each other and with the electrical signal to the source, one can determine the 

amplitude, phase and phase gradient of the field along the line of measurement. The 

locations where one of these three parameters becomes inaccurate are usually those where 
the other two parameters can be measured with enhanced precision. In this way the 

magnitude and phase of the reflection coefficient can be obtained reasonably accurately down 
to 30 Hz. 


FIGURE 12 

This figure shows the results. Although the individual points show some scatter there 
are definite trends and several peaks, or resonances are clearly evident. For example there 
is some confidence in the peaks at around 95, 130 and 200 Hz. These seismic resonances are 
consistent with the theoretical work and measurements of Sabatier, Bass and others at the 
University of Mississippi. 

In 1989 a seismic survey team drilled one or two core samples on our exact site. It 
was discovered that our site was almost ideal from an acoustical point of view. Apart from 
the top few centimeters of grass and its roots, the ground was a layer of silt of uniform 
consistency and almost constant thickness (1.6 +/- 0.3 m) lying directly over bedrock. 

FIGURE 13 

Time-of-flight measurements along the surface are shown in this figure. These were 

made by hitting a heavy metal disk lying on the ground with a hammer, and receiving the 
signal with a geophone. The sound speed in the silt layer is calculated to be 330 m/s 
(almost the same as the speed in air) and in the rock about 2000 m/s. From the break-point 
on this curve the thickness of the layer is calculated as 1.9 +/- 0.3 m. The v = 330 m/s 

part of this plot does not pass through the origin but intersects the ordinate at about t 

3.8 ms. This time delay is related to the slow sound speed through the top few centimeters 
of soil and grass-roots, but we were not able to measure the break-point due to the soil-silt 

interface. The calculated quarter- wavelength-layer-thickness resonance for the silt layer is 
about 45 Hz. 

FIGURE 14 

The acoustic-to-seismic transfer function was measured using a Mark Products L-21A 
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geophone pushed into the ground surface and a collocated microphone 10 cm above the 
surface. The two signals were analyzed and compared using a Bruel and Kjaer Model 2032 
Dual Channel Signal Analyzer. The acoustic-to-seismic transfer function was found for various 
angles of incidence ranging from normal to about 87°. Those for normal incidence and for 
84° are shown in Figure 14. Measurements at oblique incidence show a) larger surface 
admittance, b) smoother curves, and c) an upward shift in frequency by about 10%, compared 
with the admittance spectrum for normal incidence. 

Quarter- wavelength resonances in the silt layer should lead to maxima in the acoustic- 
to-seismic admittance spectrum at roughly 45, 135 and 225 Hz, and minima at 90 and 180 Hz. 
The only apparent agreement seems to be maxima at about 50 and 135 Hz and a minimum at 
about 85 Hz. Although the results could suggest a peak around 225 Hz and a dip at a 
frequency slightly greater than 180, the measurements are inconclusive. The peaks at about 
70, 105 and 180 Hz appear to be completely unrelated to the silt layer. Some of this 
structure could be due to the thin layer of topsoil and grass roots. 

Maxima in the acoustic reflection coefficient of the surface, Figure 12, should be 
related to the minima of the surface admittance spectrum. Figure 14. The match between 
these two spectra is far from satisfactory. However, the peaks of the reflection coefficient 
at about 95 and 195 Hz are not inconsistent with the dips in the admittance spectrum at 
roughly the same frequencies and are predictable, within experimental, from the thickness of 
the layer found from the core sample or the refractive survey. The peak at about 135 Hz in 
the reflection coefficient is unrelated, but the admittance spectrum does suggest a dip at 
about this frequency. 

Clearly, although our current understanding allows us to explain many aspects of these 
measurements, there are other features of this rather simple ground structure that require 
additional elucidation. Certainly more work is required before we can accurately predict the 
acoustical behavior of more realistic and complex ground structures. 
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SUMMARY 

In this paper we compare the results of Donato's exponentially varying ground model, 
Attenborough's exponentially varying ground model and the rigid backed thin layer model. We 
show that these models produce similar results for slow variations. For rapid variations the results 
are quite different but the basic theory used is only correct for the thin layer model. These results 
suggest that the exponentially varying models are not necessary for fitting ground impedance data. 

INTRODUCTION 

Donato proposed an exponentially varying ground model to be used for the interpretation of 
ground impedance data.* Attenborough has demonstrated that the exponential variation chosen by 
Donato results in model grounds with increasing porosity with depth and has derived a ground 
model which has a decreasing porosity with depth.^ 

In this paper we examine the behavior of both these models in the limit of large and small 
variation and compare the results to the rigid backed layer model. 3 To facilitate this we have 
reduced the solutions to their simplest forms and have employed Attenborough's low frequency/high 
flow resistivity results for numerical comparison. 

I. GROUND MODELS 


A . Rigid Backed Layer 

A layer of porous material of thickness d overlying an acoustically rigid surface has a surface 
impedance of the form: 


Z(0) = i Z c cot (kd) (1) 

where Z c is the impedance of a seminfinite half space of the porous material and k is the complex 
wave number in the porous material. 
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B . Donato's Exponential Model 


Donato has derived a impedance model for a material whose porosity times wave number 
decreases exponentially with depth. Attenborough has demonstrated that for natural grounds this 
implies that the porosity increases exponentially with depth and the wave number decreases 
exponentially with depth. This will not commonly occur in natural ground surfaces but may be a 
useful model in some cases. With the notation above Donato's formula becomes 


Z(0) = iZ c 


Jo(2k/a) 
Ji(2k/a) ’ 


( 2 ) 


a is the exponential varation of the square of the complex wave number 

k(z) 2 = k(0) 2 e- az . 


(3) 


C. Attenborough's Exponential Model 

Attenborough's solution for a porous material whose porosity decreases exponentially with 
depth and wave number increases exponentially with depth is given by 


Z(0) = iZ c 


Hj > 2) (2k/a) 
Hj 2) (2k/a) ’ 


(4) 


where 


k(z) 2 - k(0) 2 e az . (5) 

II. BEHAVIOR OF THE IMPEDANCE AND WAVE NUMBER 

It will be useful in the interpretation of these models to have a specific formulae for the wave 
number and impedance of a homogeneous porous material. For this paper we will use 
Attenborough's low frequency approximation: 

Z C = J&- =.218(^-Y /2 (1 + i). (6) 

yQco Vf/ 

o e is the effective flow resistivity of the material, y is the ratio of specific heats and c is the speed of 
sound in air. 2 
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m. BEHAVIOR OF THE GROUND MODELS IN THE LIMIT OF LARGE AND 

SMALL ARGUMENTS 


A . Rigid Backed Layer 

i) Limit as d — » 0. 

For a thin layer d — » 0 and Eq. (1) becomes 

(7) 


( 8 ) 

where is co/c. Note that the imaginary term approaches infinity as k Q d goes to zero, while the 
real part depends only on the layer thickness and the surface flow resistance. This form is 
displayed by Attenborough. 2 

ii) Limit as d — > oo 

As d — » °° the model should recover the result for the homogeneous half space. The 
cotangent can be expanded in terms of the exponents of the real and imaginary parts of kd. 

lim cot(kd) = lim = e:*' d e+M / eijElLejM.,- s-*» d e*** (9) 

d— d— 2 2i 


Z(0)= lim iZcCot(kd) = i^ -i^ 
d— >0 kd 3 

If we use Eq. (6) to relate Z c and k for low frequency we find 

Z (0)- 4,I (- 218 ,) 27Qdff %i-J— 

30 -ySlkod 


where 


k = kj + i k 2 . 
k 2 must be positive so that 


Z(0) = iZ c (-i) = Z c , 


GO) 


and the original condition is recovered. 
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B . Donato's Exponential Model 

i) Limit as a becomes small 


As a becomes small the medium approaches a homogeneous media. If we take the limit of 
Eq. (2) for small a and large 2k/a we find 


Z(0) - i Z c 


.J 2k<x cos /2k . K ^ 

V 2k V a 4 / 

J 27ta rn o /2k. 7t _7t \ 
V 2k ° \a 2 4) 


= i Z c cot 


/ 2k 7L \ 
Va 4 / 


(ID 


This is like the impedance of a thin layer of thickness 2/a with an additional -it/4 phase 
change. The next correction term is of order a/2k. A pressure release backed thin 
layer would have a phase change of -ic/2. As a — > 0, the cotangent term will approach -i as in 
Section A-ii) and Z(0) = Z c as expected. 


ii) Limit as a — » 

In the limit as a — > the argument becomes small and the ascending series may be used to 
evaluate the Bessel functions. 


zm “ iz <>i:- i -te • <12) 

a 

The behavior of this solution is very similar to Eq. 7. We have a rapidly increasing imaginary 
part and a constant real part as the frequency decreases for fixed d and o e . The imaginary parts are 
identical if the rigid backed layer has a thickness 1/a, while the real parts are equal if the rigid 
backed layer thickness is given by 1.5/a. 

C. Attenborough's Exponential Model 
i) Limit as a — » 0 

The asymptotic expansions can be employed for the Hankel functions giving 

. / 2 a q — i( 2 kJa - 71/4) 

Z(0) - i z c Yrc& = i z c e- i7t/2 = z c . (13) 

./2q e -i(2k/a - nU) 

» it 2k 62 6 
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The Attenborough model recovers the homogeneous half space surface impedance as a — » 0. 
ii) Limit as a — » CO 


The small argument formulae for the Hankel functions are inserted in Eq. (4) to give 


Z(0) = Z c [(|-ie)^-i^ln(|)‘, 


(14) 


where e = .5772. 


This result is not easily interpreted in terms of a layered model. The behavior of this solution 
is best illustrated by use of Eq. (6) to yield 


Z(0) = 5.923 In 


+ 3.419 + i 13.955 


(15) 


As a — » <x> the impedance of the Attenborough model has a large negative real part tending to - «> 
and a constant imaginary part This puzzling result indicates that the surface is not absorbing 
energy and has a reflection coefficient greater than one! In a gross sense the behavior is physical. 
The reflection coefficient approaches one as the impedance becomes infinite. The only problem is 
that the surface cannot be generating acoustic energy. 


iii) Limit for 2k/a > 1, a not infinite 

A third limit is developed by Attenborough as useful for computation and comparison with 
data. This form is developed for a small enough that the leading term in asymptotic series for the 
Hankel functions may be used. For 2k/a > 1 


H< 2) (2k/a) , Q+jf) . . r, , B] 

H ( , 2) (2k/a) (l ' 8 ^ 2k) 4k '* 


(16) 


and 


Z < 0) * Z c{ 1+ tk} 


Using Eq. (6) to relate k and Z c gives us 


Z(0) = Z c + 


ic Ta 

4 Y®La ‘ 


(17) 


(18) 
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Defining a e = a/Q and inserting numerical values from Eq. (6) gives us Attenborough's form: 

Z(0) = .218 (y) 1/2 + i .21 8 (y - ) 1/2 + 9 - 74 (y)] 0 9 ) 

The next terms in the asymptotic series are on the order of 7% of the last term in Eq. (19) when the 
argument of the Hankel function is one. 

Note that we can recover Eq. (13) by letting a approach zero. Also note that the second term 
in Eq. (17) is very similar to the form for the imaginary part of the impedance of a thin rigid backed 
layer. Compare 

iZ c -SL and %■ . (20) 

c 4k kd 

The second term in Eq. (17) is the imaginary part of the impedance of a thin layer of effective 
thickness d e = 4/a. The imaginary parts dominate the impedance for large a. 

IV. NUMERICAL RESULTS 

To calculate numerical values for the three impedance models we set 


kd = 2k/a = x(l + i). 

Then, using Eq. (6), we solve for f and Z c in terms of x: 

p __ axe l 2 1 
[47cyQ(.218)J ’ 

and 


( 21 ) 


( 22 ) 


47ryD(.218) 2 o e 


(1 + i) 


We use the following typical values of y, Q, o e , and a based on our experience and that of 
Attenborough: 
y = 1.4 
Q = 0.4 

a e = 120,000 MKS rayls 
a = 40. m' 1 ; d = 5 cm. 

Then, we calculate impedances using Eqs. (1), (2) and (4) for x = 0 to 5. The results are plotted in 
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Fig. 1 (rigid backed layer), Fig. 2 (Donato's formula), and Fig. 3 (Attenborough's solution). The 
imaginary parts of the impedance are multiplied by -1 so the plots of the real part are usually on the 
positive side of the vertical axis and the imaginary parts are on the negative side. The plots are 
nearly identical for values of x greater than one. For the variables above, x = 1.0 corresponds to 
654 Hz. 

Figure 4 displays the normal reflection coefficient calculated from Eqs. (1,2, and 4). The 
behavior is similar for all the models. Better agreement can be achieved between any two models 
by the choice of the equivalent depth of the exponential variation. 

V. DISCUSSIONS AND CONCLUSIONS 

The surface impedance predicted by each of the three models above approaches the 
homogeneous half-space impedance as the variation of wave number becomes small or the layer 
depth becomes large in the rigid backed model. 

As the exponential variations become larger the impedance formula can be approximated as a 
constant or slowly varying real and imaginary part plus an imaginary term which is proportional to 
cx/0) or 1/oxl. 

For very rapid variations, the expansion of Attenborough’s solution results in a non-physical 
solution (Eq. 13). 

The basic assumption in the derivation of Eq. (6) and it's more exact analogues, is that the 
gradients of the variables with respect to the propagation direction are much smaller than gradients 
of the variables normal to the direction of propagation.^ The result that the reflection coefficient is 
greater than one for small variable x is probably due to the error in Eq. (6) rather than any physical 
error in the theory leading to Eq. (4). 

By the same reasoning, Donato's formula should be inaccurate for small values of the variable 
x. There is no physical problem with the thin rigid backed layer since the porous layer is 
homogeneous and Eq. (6) should hold. For the variables we have chosen, there appears to be little 
practical reason to employ the exponential models to fit ground data, while there appears to be a 
significant theoretical reason for not using the exponential models in the region where they vary 
significantly from the rigid backed layer. 

At very low frequencies, the impedance translation theorem can be employed to calculate the 
impedance of an impedance backed layer. This model has sufficient flexibility to fit most data 
without the theoretical difficulties of the Donato or Attenborough models. 
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Figure 1. The real and imaginary part of the grbund impedance versus the parameter x for the 
thin rigid backed layer. The imaginary part is multiplied by negative one for 
display purposes. 
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Figure 2. The real and imaginary parts of the ground impedance versus the parameter x for 
Donato* s exponentially varying model. 
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Figure 3. The real and imaginary parts of the ground impedance versus the parameter x for 
Attenborough’s exponentially varying model. 



Figure 4. The absolute value of the reflection factor for normal reflection for the three 
models versus the parameter x — thin layer, Donato, Attenborough.) 
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ABSTRACT 

Measurements of acoustic pulse propagation in the 5- to 500-Hz frequency band were conducted under 
various snow cover conditions during the 1989-1990 winter in New Hampshire. The objective was to 
determine the effect of snow cover thickness and other snow properties on the absorption of acoustic pulses. 
Blank pistol shots were used as the source of the acoustic waves, and geophones and microphones in an 80- 
m-long linear array served as receivers. Snow thicknesses ranged from 0.05 to 0.35 m, and densities varied 
from 100 to 350 kg m -3 during the 10 separate measurement days. Preliminary analysis indicates that the peak 
pulse amplitude decayed in proportion to ~r -u for most conditions and that the acoustic-to-seismic ratios 
varied from about 4 to 15 x 10 -6 m s _1 Pa -1 . Theoretical waveforms were calculated for propagation in a 
homogeneous atmosphere using Attenborough’s model of ground impedance. An automatic fitting procedure 
for the normalized experimental and theoretical waveforms was used to determine the effective flow resistivity 
of the snow covers, and gave values of 10 to 35 kN s m" 4 , in agreement with earlier results. 


INTRODUCTION 

Absorption of sound energy by the ground is important in understanding noise propagation through the 
atmosphere. It affects predictions of traffic, industrial, or blasting noise levels, which are becoming 
increasingly important in mitigating or preventing community noise problems and assessing environmental 
impacts of various activities. In previous work it has been shown that a snow cover has a large effect on acoustic 
pulse propagation, causing increased attenuation and marked waveform changes compared with propagation 
over grassland (Ref. 1 ). Those measurements were for a single snow cover, so measurements were undertaken 
during the 1989-1990 winter to investigate additional snow covers and to examine the effect of snow cover 
thickness and other snow properties on pulse propagation. This paper reports on the experimental approach, 
preliminary results of data analysis, and first steps towards an automatic inversion procedure to determine 
acoustically the properties of the snow cover. 


EXPERIMENTAL MEASUREMENTS 

As in previous measurements, a .45 caliber blank pistol held and fired 1 m above the surface was used 
as the source of the acoustic waves. The receivers were a linear array of 4.5-Hz Mark Products Model L-15B 
geophones and Globe Model 100C low frequency microphones. Two Bruel & Kjaer Type 41 65 microphones 
were used to record the source pulse. Both types of microphones have a flat response in the frequency band 
of interest. A Bison Model 9048 recording system was used to acquire 48 channels of data at a 5-kHz rate. The 
bandwidth of the measurements is estimated as 5-500 Hz and is limited mainly by the source output. 
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In the fall of 1989, vertical and horizontal component geophones were installed along a relatively flat 
80-m-long line. A few geophones were also buried 0.5 m deep in the soil 30 and 60 m away from the location 
of the source. During Jhe, winter, just before each measurement period, geophones and microphones were 
installed at the snow surface and probe microphones (Ref. 2) were inserted into the snow. A number of pistol 
shot responses were then recorded, and these sensors were removed after that day’s measurements were 
completed. Only the surface sensors will be discussed in this paper. 

On the days that acoustic experiments were conducted, a snow characterization pit was dug and the 
temperature, density, grain size, and crystal type were determined for each layer present. Snow and frost depths 
were also recorded. 

Meteorological data were collected using a Campbell Scientific Model 21X data logger. Temperatures 
were measured within the ground and snow and at heights of 0.01 , 0.03, 0. 1 , 0.3, 1 , 3, and 5 m in the air. Wind 
speeds at 1- and 3-m heights were also recorded, along with relative humidity and barometric pressure. 
Measurements were taken every minute, but averages, variations (minimum, maximum, and standard 
deviation), and instantaneous values were recorded every 10 and every 30 minutes during the acoustic 
experiments. Values were recorded every 4 hours during the rest of the winter. 


ACOUSTIC WAVEFORM ANALYSIS 

Figure 1 shows the waveforms recorded on nine separate days by the Globe 1 00C low frequency surface 
microphones a distance of 60 m from the source. The positive peak amplitudes of these pulses, along with the 
air temperature, snow depth, and snow density (for the surface layer) are given in Table 1 . (Experiment 4 used 
a different sensor array than the rest of the experiments and has not yet been analyzed.) 


TABLE 1. MEASURED AMPLITUDES, ENVIRONMENTAL PARAMETERS, AND BEST FITTING 
WAVEFORM PARAMETERS FOR THE 1989-1990 WINTER EXPERIMENTS. 

Ampli- Change in 


Expt 

No 

Date 

(1989- 

1990) 

tude, 

(Pa, 

at 60 m) 

Air temp. 
(°C) 

Snow 

depth 

(mm) 

Snow 

density 

(kg/nr) 

Flow 
resistivity 
(kN s/m 4 ) 

fitted snow 
depth 
(mm) 

1 

29 Dec 

3.1 

-12.4 

185 

170 

25 

0 

2 

4 Jan 

4.9 

3.1 

170 

260 

30 

0 

3 

10 Jan 

4.3 

1.3 

140 

280 

35 

-50 

4 

19 Jan 

17.0 

-3.0 

50 

210 

— 

— 

5 

22 Jan 

2.0 

-5.3 

190 

100 

10 

+50 

6 

31 Jan 

2.2 

-2.8 

350 

140 

10 

0 

7 

8 Feb 

1.9 

3.0 

280 

150 

10 

+50 

8 

6 Mar 

3.3 

-4.0 

140 

340 

35 

-50 

9 

15 Mar 

16.1 

14.3 

0-60 

350 

— 

— 

10 

12 Apr 

16.7 

3.2 

0 

300 

— 

— 


Note: The snow cover was continuous for all of the experiments except for experiment number 4 (9/ 10 of the 
ground was covered), experiment number 9 (5/10 covered), and experiment number 10 (no snow). 


52 


The two largest arrivals were recorded on days when there was little or no snow cover present, and have 
amplitudes about five times larger than the pulses recorded when snow was present. The waveforms recorded 
over snow are all elongated to various degrees, and exhibit relatively stronger low frequency content than those 
recorded without snow present. 

In Reference 1, a method of calculating theoretical acoustic pulse waveforms from known surface 
properties was developed and verified. The procedure is briefly outlined here. For a monofrequency source 
in the air and a receiver on the surface, the acoustic pressure a slant distance (r) away from the source is given 
by 


P/P 0 = l/kr e* r (1+0) 

where P Q is a reference source level, k is the wave number in air, Q is the image source strength representing 
the effect of the ground, and <r'°* is suppressed. At high frequencies (kr » 1), Q can be written as (Ref. 3 
and 4) 


Q=R p + (\-R p )F(w) 

where R is the plane wave reflection coefficient, F is the ground wave term, and w is a numerical distance, 
all of which depend on the specific impedance Z, of the ground. The impedance is itself dependent upon 
frequency; thus, so is Q. [The elongation and relatively stronger low frequency content of the measured 
waveforms in Figure 1 can be explained theoretically by the decrease in R p at high frequencies and the 



Figure 1. True amplitude, time aligned, low frequency 
surface microphone waveforms at 60-m range 
from a .45 caliber pistol shot 1 m high above 
the snow or ground surface. These waveforms 
were recorded with the same microphone on 
nine separate days, and the numbers refer to 
the measurement days listed in Table 1 . The 
two largest waveforms occurred on days when 
there was very little or no snow cover present. 
Note that waveforms 3 and 8 are slightly mis- 
aligned in time; the shift is the result of the low 
frequency portion of the waveform being larger 
than the direct arrival. 
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enhancement of F(w) at low frequencies (see Ref. 1, Fig. 4).] By determining Q over the frequency band of 
interest, an inverse FFT*can be used to construct theoretical pulse waveforms in the time domain. Nicolas et 
al. (Ref. 5) have shown that an explicitly layered model of the ground must be used to represent thin snow 
covers, and this was done in the calculations presented here using 

Z = (Z 3 - / Z 7 tan k n d) / (Z 0 - i Z 3 tan k^d) 

where d is the snow layer thickness, k 0 is the wave number in the layer, and Z 7 and Z 3 are the impedances of 
the layer and substratum, respectively' (Ref. 6). 

The impedance Z 7 and wave number k^ of the snow were calculated using Attenborough’s (Ref. 7) four- 
parameter model. For all of the calculations,’ the grain shape factor n' was set to 0.5 and the pore shape factor 
ratio s ( was 0.8. The porosity Q was determined from the measured density of the snow, and the effective flow 
resistivity a was allowed to vary. 

A new result presented in this paper is a method of comparing calculated and observed acoustic pulse 
waveforms. A suite of waveforms were calculated and the best fitting waveform was selected under the L, noma 
criterion (i.e., the sum of the absolute value of the differences between the calculated and observed waveforms 
over a fixed time window). A least squares criterion, the L 1 norm, was avoided because it heavily weights, and 
tries to reduce, the maximum misfit. Since the source pulse in the calculations is an assumed one, and not ac- 
tually measured, I wanted to allow for errors in this assumed pulse to be ignored while accurately fitting the 
overall, low frequency portion of the measured waveforms accurately. 

Eight theoretical waveforms were calculated to fit the observed waveform at r = 60 m using the measured 
snow thickness and porosity, with the effective flow resistivity o varying from 5 to 40 kN s m . Then, for the 
best o, four additional waveforms were calculated, with the snow thickness changed by ±0.05- and ±0.1 -m 
increments from the measured thickness to see if the fit could be improved. An example is given in Figure 2. 

♦Fast Fourier Transform 



d + O.fm 
d -+ 0.05m 
d — O.QSrn 
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0 . 02 * 


Figure 2. Comparison between normalized measured 
and calculated waveforms for experiment 
number 6 (see Table 1 ) at a range of 60 m. The 
solid lines are the measured waveform; the 
dashed lines are calculated waveforms with 
the indicated effective flow resistivities o. 
The measured snow depth d was 0.35 m. Stars 
mark the best fitting waveform. 
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0.02s 



Figure 3. Comparison between normalized meas- 
ured (solid) and calculated (dashed) wave- 
forms for experiment number 6. The wave- 
forms at all the ranges were calculated using 
the parameters from the fitting procedure at 
60 m. At 10-m range, a Bruel & KjaerType 
4 1 65 microphone 0.3 m above the snow was 
used (and the measured waveform shows 
some evidence of being clipped); the other 
measurements were made with Globe Model 
100C low frequency microphones on the 
snow surface. 


Figure 4. Comparison of normalized measured 
(solid lines) and calculated (dashed) wave- 
forms at 60 m for propagation over various 
snow covers. The numbers refer to the ex- 
periment numbers given in Table 1 , where 
the best fitting flow resistivities and snow 
layer thicknesses are listed. 


In this case the best fit was obtained for o = 10 kN s irf 4 and the measured snow thickness. Using these best 
fit values of o and d, more waveforms were then calculated for different propagation ranges. The comparisons 
between these waveforms and observations are shown in Figure 3, and the agreement is quite good. 

All the measured and best-fit calculated waveforms for snow are shown in Figure 4. The fitting 
procedure has been able to automatically match waveforms of quite different appearance. The last two columns 
of Table 1 list the effective flow resistivities and snow depths determined using this fitting procedure. In all 
cases the snow thickness was within ±0.05 m of the measured thickness, a reasonable variation considering 
the variation in the actual snow cover thickness across the propagation path. 


ADDITIONAL ACOUSTIC MEASUREMENTS 

The amplitude decay as a function of range was determined by least squares fitting of the data from the 
low frequency microphones to 

A(r) = A(r 0 ) r- 

where r is the propagation distance in m, A(r) is the peak amplitude in Pa at range r,A( r () ) is the source amplitude 
at a reference distance r , and a is the distance attenuation exponent. For the data analyzed so far, the results 
are given inTable 2. Values of a for snow range from 1.6to 1.9, compared with the expected 1.0 from spherical 
spreading. For the last two experiments, with little or no snow present, the coefficient is around 1.1. 
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TABLE 2. RANGE DECAY COEFFICIENT AND ACOUSTIC-TO-SEISMIC COUPLING RATIO 
MEASURED FOR AIR WAVES. 


Range decay Acoustic-to-seismic 

coefficient a coupling ratio, m s _l Pa -1 


Expt 

No. 

Date 

No. of 
points 

Value 95% 
of confidence 
a interval 

No. of 
points 

Value of 95% 

coupling confidence 
ratio interval 

1 

12-29-90 

12 

-1.60 ± 0.30 

10 

4.01 ± 1.25 x 10^ 

2 

1-04-90 

13 

-1.71 ± 0.49 

10 

6.35 ± 5.18 

3 

1-10-90 

18 

-1.69 ± 0.22 

15 

3.41 ± 0.87 

4 

1-19-90 





5 

1-22-90 

8 

-1.76 ± 0.64 

6 

5.10 ± 5.03 

6 

1-31-90 

18 

-1.84 ± 0.22 

15 

15.1 ± 2.77 

7 

2-08-90 

16 

-1.91 ± 0.36 

13 

5.93 ± 1.49 

8 

3-06-90 

18 

-1.73 ± 0.16 

12 

4.27 ± 1.07 

9 

3-15-90 

11 

-1.05 ± 0.46 

9 

6.28 ± 2.66 

10 

4-12-90 

30 

-1.12 ± 0.19 

25 

10.2 ± 1.69 


The ratio of induced particle velocity in the snow or soil to incident pressure was determined from the 
collocated surface vertical component geophones and surface microphones (Table 2). These ratios vary from 
3 to 15 x 10 -6 m s -1 Pa -1 . Note that some of the values have very poor confidence intervals (e.g., experiments 
2 and 5). It is hoped that these values will be better determined when all of the data are analyzed. 


CONCLUSIONS 

The experiments were successful in obtaining accurate measurements of pulse propagation over a 
variety of seasonal snow covers. Preliminary values have been presented for the range decay and acoustic-to- 
seismic coupling coefficients, and more accurate values will be provided when the data analysis is completed. 
I have also demonstrated a waveform matching procedure that can be used to select the theoretical waveform 
that best fits the measured data. 

Future work will include completing the data analysis, including determination of attenuation coeffi- 
cients in the snow from the probe microphone recordings, and correlating the acoustic effects with the snow 
cover properties. A true waveform inversion procedure will also be developed. 
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SUMMARY 


A calculation method is presented for sound propagation over an impedance 
discontinuity in flat ground with a homogeneous, still atmosphere. The method is 
based on an approximate solution to a two dimensional boundary integral equation 
formulation of the problem, which expresses the wave field as the solution for 
homogeneous ground plus an integral over half of the boundary. Through recognising 
this integral as a generalised Fourier integral, asymptotic methods are applied to 
evaluate the part of the integral most expensive to compute by numerical quadrature. 
Single frequency excess attenuation results for propagation from a point source 
above rigid ground to a receiver above absorbing ground are discussed. The results 
are applied, with air attenuation and A-weighting, to a notional jet engine noise 
source; simple trends are noted. 


INTRODUCTION 


The problem discussed in this paper is propagation from a point source in a 
homogeneous still atmosphere above flat locally reacting ground. Efficient 
calculation methods for the wave field above acoustically homogeneous ground are 
well known (e.g. ref. 1). More recently sound propagation over impedance 
inhomogeneities has been theoretically examined; a thorough review is given in 
reference 2. A limitation of the accurate calculation methods is their 
computational expense . 

Here we focus on propagation over a single straight line impedance discontinuity 
which lies perpendicular to the direct source-receiver propagation path. A 
development to an existing calculation method is described which significantly 
reduces the computational expense. 

The improved calculation method is applied to grazing incidence propagation from 
a source above a rigid surface to a distant receiver above absorbing ground. 
Monofrequency excess attenuation results are examined and some simple trends are 
observed. The results for a l-5m high receiver are applied, with air attenuation 
and A-weighting, to a notional jet engine noise source at l*5m height. Again some 
simple trends are noted. 
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CALCULATION METHOD 


t 


Description of the Problem 

Figure 1 illustrates the problem. A point source with harmonic time dependence 
( e ~io)ty i s situated over a flat locally reacting surface of infinite extent. The 
surface is divided by a straight line into two half planes. Each half plane is 
acoustically homogeneous and characterised by a frequency dependent complex 
admittance (the inverse of the normalised acoustic surface impedance). We are 
interested in evaluating the acoustic potential at a point in a vertical half plane 
that is bounded by the surface, passes through the source and is perpendicular to 
the line of the admittance discontinuity. For the mathematical description we will 
use right-handed Cartesian coordinates Oxyz as indicated in Figure 1, the y-axis 
vertical and the surface in the plane y=0 . The source and receiver coordinates are 
(0,h s ,0) and (L,h r ,0) respectively. The admittance discontinuity is along the line 
x=X in the surface . 


An Existing Calculation Method 


First we consider the related problem in which the source is replaced by an 
infinitely long coherent line source, parallel to the admittance discontinuity. 
This cylindrical wave propagation problem is mathematically equivalent to the two 
dimensional problem which is illustrated in Figure 2. From the mathematical 
expression of this problem as a two dimensional boundary value problem (consisting 
of the Helmholtz equation and suitable boundary conditions) the following boundary 
integral equation can be derived (ref. 3): 


- G ( 32 (t 1 ,t 2 ) + ik((5,-p 2 )J* oo *(s, t, )G 0 2 (s, t 2 )dx. 


( 1 ) 


In this equation ^=*(0,1^) is the source position, t 2 =(L,h r ) is the receiver 
position, and /3 1 and (3 2 are the admittances of the two halves of the boundary. The 
integration is over the interval ^=(-00 ,X] ; this is the part of the boundary with 
admittance /3 1 . s=(x,0) is a point in the boundary. For two points a and b, <t>(a,b) 
denotes the acoustic potential detected by a receiver at a when insonified by a unit 
source at b; G^(a,b) (where |3= = /3 1 or (3 2 ) denotes the same quantity in the simple 
case when the boundary has homogeneous admittance j3. Efficient methods for 
evaluating the solution in this simpler case have already been developed (refs. 4,5). 


Equation (1), which describes an inhomogeneous admittance boundary problem, can be 
solved accurately for by the boundary element method (refs. 3,4,6). We 

consider here an approximate but less computationally expensive method of solution. 
To develop this we make the physically plausible assumption that the potential in y 
is what it would be if the whole boundary had admittance |3 1 (refs. 7,8). Thus 
$(§,£i) equation (1) is replaced by G^^s,^), giving the following approximation 

to <fc(t 2 , t, ) : 


where 


4 > a(£ 2 ,£i) = G 02 <£i ,t 2 ) 


+ i(0,-0 2 )I(X) 


•x 


I(X)= k 


G /S i < S . t n 2 ( § , t 2 ) dx • 

—00 


( 2 ) 


This approximation avoids using the boundary element method. 


60 


Using < 1 , a(£2>£i) we can calculate Qa(£ 2 >£,)> an approximate cylindrical wave 
reflection coefficient. Let d and D denote the distances from source and from image 
source (at (0,-h s )) to receiver, respectively. Assuming a source with unit volume 
flow rate amplitude, 

*a(£ 2 ,£i> - -^HjO(kd) + Q A (£ 2 ,£ 1 )[-jH^)(k D )] , 

where H^ 1 ) is the Hankel function of the first kind of order zero. 

For propagation over short distances, the difference 201og 10 itf’A^ ,t )l 
-201og^Ql < * ) (£ 2 .£ 1 ) I . where ^(tj.t^ is calculated by the boundary element 1 method, has 
been found to be around 0-ldB (ref. 6). This suggests that is an accurate 

approximation to Q(t 2 ,t,), the exact cylindrical wave reflection coefficient for the 
two dimensional problem illustrated in Figure 2. Also, it has been argued that 
(ref. 6), for a receiver in the far field of the image source (kD>l) , Q(t 2 ,t ) is an 
accurate approximation to the spherical wave reflection coefficent, q, for the three 
dimensional problem illustrated in Figure 1. If the point source in the three 
dimensional problem has unit volume flow rate amplitude, then the acoustic potential 
at the receiver position is 

e ikd e ikD 
^ 4ird ^ 4irD 

Replacing q with Q(t 2 ,t,), which is approximated by Qa(£ 2 .£,)> we obtain an 
approximation for <p : 

e ikd e ikD 

*A = - -2£d " Q a(£2,£,)-4^d • (3) 

We assume throughout the rest of the paper, without further comment, that this is a 
good approximation for <p. 

The main computational expense in this approximate calculation method is in 
evaluating the integral I(X). In previous calculations (ref. 6) I(X) was evaluated 
numerically after first replacing the lower limit of integration, - 00 , by a 

^• c ^- an ^'3.y large negative value. Unfortunately the integrand in equation (2) is 
usually highly oscillatory over the range of integration, making numerical 
integration an expensive process. Here we derive a semi-analytical method of 
evaluation which deals efficiently with the part of the integral that is most 
expensive to evaluate numerically. 


The Improved Calculation Method 

We begin by examining the general behaviour of the integrand in equation (2). 

We note that, for a receiver at the point s in the boundary, we can write 

4iG /S(§,t) - 1 ) (k| t-s | )R|j(s , t) , 

where Rp(s , t)=l+Q^(s , t) , and Q^(s,t) is the cylindrical wave reflection coefficient 
for a homogeneous surface of admittance 0. As x is increased from - 00 , the real and 
imaginary parts of H^O(k|t-s|) oscillate in a well defined fashion, while Ro(s,t) 
changes less rapidly. In fact, from the asymptotic expansion of the Hankel function 
at large argument, we see that if we factorise, 

G /3 ( s,t) = eik|t-si s(§i t) , 
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then S(s,t) is smooth as a function of x compared to ei k| £ - § ( when k|t-s| is large. 
This observation suggests that I(X) can be usefully written in the form of a 
generalised Fourier integral: 


I(X) =[ X f(x)e ik S( x )dx , 
** —00 


( 4 ) 


f(x) = kGp, (s,t,)G|5 2 (s,t 2 )e i k g( x ) , 
g (x) = g, (x) + g 2 (x) , 

g,(x) = |t,-S|, g 2 (x) = lt 2 -s|. 

Notice that g(x) is the distance from source to receiver via the point s. The 
location on the boundary of the geometrical reflection point, x r , is therefore given 
by g'(x r )= 0 . When s is sufficiently distant from t 1 , t 2 , and (x r , 0 ), f(x) is a 
slowly changing function of x compared to eikg(x) . This fact allows us to use 
simple asymptotic methods to help evaluate I(X). 


where 


and 


To introduce the asymptotic analysis we consider first what proves to be the 
simplest type of configuration to deal with. This has x r well outside 7, and and 
t 2 at least one wavelength from the boundary. For this type of configuration we may 
integrate I(X) by parts to give 

I(X) = J,(X) + R,(X) ( 5 ) 

where 


J,(X) = fm e lk g( x ) 


ikg'(X) 


R,(X) 


[ x d_ 

f(x) ^ 

dx 

—00 

ikg' (x)J 


: ikg(x) dx 


and then integrate R, (X) by parts to give 

R,(X)-J 2 (X)+R 2 (X) . 

where 


J 2 (X) = 


g"(X) f'(x)l f(X)e ik e( x ) 

g’lxy f (X) J ( ikg’ (X) ) 2 


R 2 (X) 


f X d_ 

d 

f (x) 1 

1 

dx 

—00 

dx 

ikg’ (x) 

ikg' (x) 


e ikg(x) dx 


If f(x) were completely independent of k, we could apply the Riemann Lebesgue 
lemma to show that, for n~l , 2 , 

R n (X) = o(k~ n ) , k -> 00 . 

In fact f(x) depends weakly on k, but f(x) approaches a limit independent of k as 
k->oo with other variables fixed (ref. 2 , p. 585 ). Thus I(X) has the following 
asymptotic approximation: 

I (X) - J , (X) + J 2 (X) , k ^ 00 . 
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J 1 (X) and J 2 (X) are the first and second terms in what is approximately an 
asymptotic expansion of I(X) in inverse powers of k. When k is large enough, 
J 1 (X)^J 2 (X), so that we can safely approximate 

I (X) * J,(X). (6) 

The above arguments do not tell us how large k should be in any particular case 
for approximation (6) to be valid. However, it is plausible that R 2 (X)<J 2 (X) when 
^ 2 (X) i (X) . Thus we can estimate the relative error in approximation (6) by the 
following upper bound on i J 2 (X)/J 1 (X) | : 

E r (X) 


f g" (X) 3 

r_L_ + 1 1 

1 1 

llg’ (X) I 21 

Lg, (X) g 2 (X)J 

Jklg' (X) 1 


(To obtain this expression, |f*(X)/f(X)| has been replaced by 

(3/2) [l/g 1 (X) + l/g 2 (X)], which is expected to be an upper bound on |f f (X)/f(X)| in 
all cases (ref. 2, p.598).) We can estimate the absolute error in approximation (6) 
by the following upper bound on | J 2 (X) | : 

lf(X)| 


E a (X) - E r (X) 


klg f (X) | 


Both E r (X) and E a (X) are infinite at X=x r and tend to zero as X-»-oo. Moreover, a 
graphical examination of E r (x) and E a (x) suggests that they are monotonic in 
(-°o,x r ) , for typical geometries, admittance values and frequencies. 


We move on to consider configurations for which still X<x r but x r -X is small 
enough for E r (X) and/or E a (X) to be unacceptably large. For the moment we require 
that both the source and receiver are many wavelengths above the boundary. The 
following breakdown of the integral is used: 

I (X) - I (t ) + K , (7) 


where K is the integral over a truncated interval 7T=U>X], 

f X 

K - kj ^/? 1 (§>J 1 )G|3 2 (s,t 2 )dx , 


which will be evaluated numerically. We will choose r so that we can satisfactorily 
approximate 

KO - J,(t) . (8) 


To reduce the expense in evaluating K numerically we want to choose r as close 
to x r as possible while still insisting that approximation (8) should satisfy 
certain relative and absolute error criteria. We can uniquely define two upper 
limits, t r and r a , for r by 

E r (r O — 0 * 2 , 


E a (r a ) - e , (9) 

where e is an arbitrary positive constant. r<r r ensures that J 2 (r) is sufficiently 
small compared to J n (0 for E a (r) to be an accurate estimate of the absolute error 
in approximation (8). Therefore, if also 7 < 7 a , then the absolute error in the 
approximation (8) is <e . Thus by taking r to be the minimum of r r and r a we ensure 
that the error made in replacing I(r) by J^r) in equation (7) is <c . 
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We can now consider the more general configuration in which x r may be anywhere 

in relation to X, but the source and receiver remain many wavelengths above the 

boundary. If X is less than x r , one of the above calculation methods applies. If X 
is greater than x r , reciprocity can be invoked (reflect the problem in the plane 
x=L/2, then swap the source and receiver) and then one of the above methods applied. 

Finally we note why the source and the receiver have so far been kept at least 

one wavelength above the boundary. If t 1 (or t 2 ) is very close to y the 

approximation (6), which involves neglecting the integral R 1 (X) in equation (5), 
breaks down. This is because, for x in a small range of y around t 1 (or t 2 ) , f(x) 
changes rapidly with x. Thus the derivative of f(x) in the integrand of R, (X) is 
very large . 

To avoid the consequent inaccuracies which may occur when t n or t 2 is within one 
wavelength of y, an additional criterion is used for the choice of r . We require r 
to be small enough so that the line x<r in y is always at least one wavelength from 
1 1 and 1 2 . 


GRAZING INCIDENCE RESULTS 


Monofrequency Excess Attenuation 

We can use the method described above for estimating ^ to examine propagation 
over flat ground through a homogeneous still atmosphere. The monofrequency excess 
attenuation over geometrical spreading due to the presence of the ground can be 
approximated by 

ikd 

Al - 201og, 0 | 5 S a -^-| d» . (10) 

Propagation from a source above rigid ground to a receiver above absorbing 
ground has been examined. To model this problem |3 1 was set to zero and the 
dependence of (3 2 on frequency was calculated by the Delany and Bazley semi-empirical 
formula (refs. 9,10), with an effective flow resistivity of 10^kgs -1 m” 3 . This value 
was chosen as being a low value for grassland (ref. 11) . (It is found that using 
two or three times this flow resistivity value causes only a small reduction in the 
magnitude of the A^ results, and no change in the trends was observed.) Six 
configurations of source and receiver heights (h s and h r ) were examined. 

Specifically, heights of 5m, l*5m, and 0-5m were used, with h s >h r in all cases. For 
each h s and h r combination calculations were carried out at four distances: L=250m, 
500m, 1km, and 2km. 

We will examine the significance of the proportion of rigid ground between the 
source and receiver. We can define a useful variable, p r , by 

p r = X/L . 

When the impedance discontinuity is between the source and receiver (0<X<L) , p r 
gives the proportion of rigid ground between the source and receiver. 

A sample of the results examined is shown in Figures 3a and 3b. Plots like 
those shown were calculated for all the octave band centre frequencies between 100Hz 
and 5kHz. Notice that, as is of course expected intuitively, when p^XO or p r >l the 
modelled ground behaves as an acoustically homogeneous plan, absorbing or rigid, 
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respectively. We therefore now concentrate the investigation only on the ranee 
0<p r <l. 6 

Figure 4(a) illustrates a simple curve shape that occurs whenever both kh s and 
kh r are small enough. Half of the plots examined were of this type. We see that Ai 
increases monotonically with p r . The gradient of each curve is greatest at the ends 
of the range of interest, i.e. near p r =0 and p r =l. At low enough frequencies the 
curves straighten out. 

Figure 4(b) illustrates the disintegration of the orderly patterns seen in 
Figure 4(a) that occurs when h s , h r , or the frequency is increased sufficiently. 
About a third of the plots examined showed this type of disorderly pattern. 

Figure 4(c) illustrates a different pattern that sometimes occurs when the 
source is higher than the receiver, but neither are so high above the boundary that 
the disorderly pattern seen in Figure 4(b) occurs. In Figure 4(c), A 1 is less 
dependent on the location of the impedance discontinuity when p r <0-3. Notice that 
the right half of this plot shows the features observed in Figure 4(a). About a 
sixth of the plots examined showed this pattern. 

When the orderly patterns seen in Figure 4(a) and the right hand side of Figure 
4(c) occur, there is usually a range of octave band centre frequencies around 1kHz 
at which some or all of the curves on a plot are separated, in most of the range 
0<p r <l, by about 3dB. This approximate 3dB increase in A^ per doubling of L when 
propagation is over an admittance discontinuity occurs only for a range of values of 
kjL- The start of this inhomogeneous absorbing ground effect corresponds with the 
start of a 6dB separation of the A x curves at p r <0, which occurs when kL=2300. This 
6dB increase in A]_ per doubling of L is a homogeneous absorbing ground effect which 
has been predicted theoretically (ref. 12). Figure 5 illustrates these 
observations . The range of plots examined show that the inhomogeneous ground effect 
fails to occur when kh s or kh r is large. This failure is observed in the left hand 
side of Figure 4(c) . 


Jet Engine Noise 

We move on to examine the excess attenuation of a notional broad band 
environmental noise source. A simple spectral shape representative of a jet engine 
Full thrust is chosen. The free field lm third octave band sound pressure level 
is taken as constant up to 200Hz, above which frequency it is reduced by 0-8dB per 
third octave band. Third octave band excess attenuations due to the presence of the 
ground are approximated here by A]_ values given by equation (10), using the band 
centre frequencies. To make the calculation more realistic we include the B.S.5727 
(1979) third octave band free field air attenuations for 20 *C and 70% relative 
humidity . 

We consider only source and receiver heights of l-5m so that the simple A^ 
pattern illustrated in Figure 4(a) dominates the results. The process of intensity 
addition over the third octave bands will produce more moderate excess attenuations 
for the broad band noise than those calculated for monofrequency sound. 

The excess attenuation of the broad band source noise caused by the presence of 
the ground and by air absorption along the propagation path is 

A 2 “ Sx - (S 2 + 201og, 0 L) dB(A) , 
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where is the total A-weighted sound pressure level of the notional source at lm 
distance in the free field, and S 2 is the total A-weighted sound pressure level that 
we calculate at the receiver position. 

Figure 6 shows how A 2 depends on p r and L. We can see that it is not possible 
to predict the A 2 values at intermediate p r from a linear interpolation between the 
values at p r =0 and p r =l . Notice however that, in the range 0 * 25<p r <0 * 75 , A 2 appears 
to vary linearly with p r . Also, in this range, A 2 increases by 5dB(A) per doubling 
of L. This dependence on L is in reasonable agreement with a commonly used 4dB(A) 
extra attenuation (of perceived noise level) per doubling of receiver distance from 
an aircraft source very close to the ground (of unspecified admittance) (ref. 13). 

At L=lkm and 2km, A 2 varies linearly with p r in the range 0<p r <0-75. This is 
useful because it means that A 2 (p r ) can be estimated from A 2 (0), the value for 
homogeneous absorbing ground, which is easier to calculate. The simple predictive 
equation, which is shown in Figure 5 for L=lkm and 2km, is 

A 2 (p r ) “ A 2 (0) - 14 • 7p r dB(A) , 0 < p r < 0-75, (11) 

Unfortunately we know of no practical results with which to compare this equation. 


CONCLUSIONS 


An improved calculation method has been presented for sound propagation over a 
straight line impedance discontinuity in flat ground. The method is restricted to 
the case when the impedance discontinuity is perpendicular to the direct source to 
receiver propagation path. The method is derived from an asymptotic analysis at 
large wavenumber of an approximate solution of a two dimensional boundary integral 
equation. Accuracy is adequate for the purpose of examining environmental noise 
propagation in ideal conditions. A limitation of the method is the assumption of 
homogeneous still air and flat ground. 

Results for long distance grazing incidence monofrequency propagation show that 
the dimensionless heights (height multiplied by wavenumber) of the source and 
receiver above the ground are as important as the location of the impedance 
discontinuity. When these dimensionless heights are small enough, the results are 
very orderly, as illustrated by Figures 4(a) and (c) . 

Theoretical results for the excess attenuation, including air absorption, of a 
broad-band A-weighted notional environmental noise source have been examined. A few 
simple trends have been noted, in particular, equation (11). 
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Fig. 4a Monofrequency excess attenuation over geometrical spreading at four long 
distances, plotted against the proportion of rigid ground. h s “h r «l*5m, f“250Hz. 
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Fig. 4b Monofrequency excess attenuation over geometrical spreading at four long 
distances, plotted against the proportion of rigid ground. h s =5m, h =1.5m, f=2kHz 
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INTRODUCTION 

An acoustic detection range prediction model (ADRPM-VII) has 
been written for IBM PC/AT machines running on the MS-DOS operating 
system. The software allows the user to predict detection distances 
of ground combat vehicles and their associated targets when they are 
Involved in quasi-military settings. The program can also calculate 
Individual attenuation losses due to spherical spreading, 
atmospheric absorption, ground reflection and atmospheric refraction 
due to temperature and wind gradients while varying parameters 
effecting the source-receiver problem. The purpose of this paper is 
to examine the strengths and limitations of ADRPM-VII by modeling 
the losses due to atmospheric refraction and ground absorption, 
commonly known as excess attenuation, when applied to the long range 
detection problem for distances greater than 3 kilometers. 


BASIC ASSUMPTIONS OF ADRPM-VII 

The basic assumptions of ADRPM-VII are the following: 

o ADRPM is based on simplified atmospheric conditions adjusted 
to a standard day during the seasonal year. In the real world, a 
standard day does not exist since temporal variations must be 
allowed for in all environmental propagation measurements. The 
effect of these variations can only be measured with sound speed 
profile soundings. 

o The noise emitted by the source is omnidirectional, broadband 
and continuous. 

o The primary propagation path is near the surface of the 
ground . 

o All attenuation elements are considered independent of each 
other with the total attenuation arrived from the summation of its 
individual parts. 

o The ground is defined as a rigid plane or a plane of finite 
impedance and the model uses a table of values of ground cover loss 
that is linearly dependent on the distance from the source. 

o " The model is developed in the context of a need to estimate 
noise levels of surface vehicles at distances ranging from tens of 
meters to hundreds of meters for a relatively wide range of 
environmental conditions” according to Fidell and Bishop (ref. 1). 
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ATTENUATION DUE TO REFRACTING ATMOSPHERES 


The model calculates propagation loss in a refractive 
atmosphere by applying a correction term to the reflected and 
surface wave terms derived from non^ref racting atmospheres. This 
correction term, which is based on ray tracing, considers the 
existence of shadow zones for upward refraction and an intensity 
ratio modification for the downward refracting case (ref. 2). 

Several representative atmospheres have been chosen from the 
given meteorological profiles in ADRPM for analysis of the models 
refractive effects. Average wind velocities u(r), surface roughness 
parameter z(o), and Monin stability length L are given for each 
selected profile: 

Neutral Profiles : Vertical temperature lapse of -.01 degrees Kelvin 

per meter and turbulence due to wind only. The following latitude 
and season was chosen for analysis: 

1. Mid-latitude (45°N), summer, with 

u(r) = 3.3 mph, 
z ( o ) = 0.15 

surface temperature = 73.8°F. 

Stable Profiles : A positive temperature gradient and damped 

turbulence due to thermal inversion only. 

1. Mid-latitude (45°N), summer night, with 

u( r) = 2.5 mph 
z(o) = 0.15 
L =39.65 

surface temperature = 62°F 
temperature gradient = .02 

= .01 

2. Midlatitude (45°N), winter night, with 

u(r) = 4.4 mph 
z(o) = 0.15 
L = 38.6 

surface temperature = 21°F 

temperature gradient = 0.07 for 0-40 meters 

= 0.02 above 40 meters 

Unstable Profiles: 

1. Midlatitude (4 5°N) , summer daytime, with 

u(r) = 3.6 mph 
z(o) = 0.15 


L = -16.88 
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2. Midlatitude (45°N), winter daytime, with 


u(r) = 6.5 mph 
z(o) = 0.15 
L = -243.5 

surface temperature = 36°F 
temperature gradient = -.02 

= -.01 

= -.004 


for 0—15 meters 
1 5—2 5 meters 
above 25 meters 


A:Attenuation Due To Upward Refraction 

The upwardly bending sound energy algorithms have evolved 
through the efforts of several investigators, with Felt (ref. 3) 
making the greatest contribution. Felt's ray tracing procedure 
requires a numerical solution to a differential equation to 
determine the ray path as a function of the initial angle of 
propagation. For a specified source height h(s) and receiver height 
h(r), attenuation is based on the distance to the shadow zone d(s), 
which is defined by: 

d(s) = ( h(s)/k ) 1 / a + ( h(r)/k ) 1 ^ a ( 1 ) 

where: h(s) = source height 

h(r) = receiver height 

d(s) = distance to the shadow zone 

and a,k are parameters that are determined from Snell’s law of 
refraction for various meteorological profiles. 

The attenuation due to upward refraction is capped by a maximum 
frequency dependent value that is dependent on the distance to the 
shadow zone, as determined from equation 1. The value of attenuation 
A(e) is calculated from: 

A(e) = A( max) ( 1- d(s)/d ) (2) 

For a source to distance receiver d, the model considers two cases: 

d < d(s) where the receiver is not in the shadow zone 

d > d(s) where the receiver is in the shadow zone 


B : Attenuation Due To Downward Refraction 

For the downwardly refracting case, a fitting function based on 
the initial propagation angle oC and the distance from the source 
to where the ray strikes the ground x is given by (ref. 4): 

tan dL = MX b (3) 

where M,b are determined in much the same way as a,k were determined 
for the upwardly refracting case in equation 1. 
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ATTENUATION DUE TO GROUND IMPEDANCE 


The attenuation due to the effect of a sound wave interacting 
with a surface of finite impedance is based on the work by Embleton, 
Piercy, Olson (ref. 5) and Delany ,Bazley (ref. 6). ADRPM-VII 
calculates the effect of ground impedance based entirely on the 
coherence of incoming waves. However, the stable conditions assumed 
for the phase dependent calculations are unlikely to exist for 
longer ranges since the effect of inhomogeneity on the delicate 
phase relationships is ignored. 

Nevertheless, the theory predicts losses of 50-70 dB for some 
conditions. Since losses beyond 30 dB are rarely observed, the model 
handles this empirical discrepancy by decreasing the effects of 
ground impedance for distances greater than 500 meters. 

In addition , the model accounts for a non-uniform surface by 
requiring a single user supplied parameter. This parameter, h, is 
the root mean square surface roughness height. Based on reference 6, 
h yields a smoothness, s, that represents the fraction of the 
reflected energy that is specularly reflected. 

However, the unique topography along the propagation path is 
not included in the model. This is an important omission since 
sloping ground can control the phase as well as serve as a barrier 
by intercepting incoming rays. 


RESULTS AND DISCUSSION 

Field data of stationary and moving helicopters have been 
analyzed over ranges from 300 meters to 12 km. The results show a 
built-in variability of the continuously received signal for ranges 
between 2 and 5 km. At these source-receiver distances, the 
refractive atmospheric state, with all its existing temperature and 
changing wind directions, will have a variable attenuation effect on 
the propagating rays and consequently produce a variable received 
signal. 

In the field, it remains difficult to determine the unique 
local sound speed profile for all threat directions, especially since 
the sound speed profile can change with the next gust of wind or the 
next reversal of wind direction. This problem of measuring time 
varying speed profiles occurs at all field locations that we have 
visited across the United States. However, the meteorological 
conditions are still determined only at the detector during ground 
vehicle testing. 

The area of the atmosphere that primarily effects ground 
vehicle vulnerability for the medium detection distances is in 
constant change due to its turbulance. A wave propagating through 
this boundary layer is variable in amplitude and is influenced by 
the daily cycle of stable and unstable meteorological conditions 
that repeat themselves several times each day. TACOM data shows that 
noon time provides the largest variation of amplitude, sometimes as 
much as 7 to 8 dB. The fluctuations are less and also slower during 
the morning and early part of the evening. In all cases, it is best 
to obtain sound speed profiles each time that a set of data is 
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measured, with as many locations as possible, but at least two 
extreme readings that would cover the source and the projected 
receiver distance* 

For ranges beyond 5 km, field data signals are intermittent, 
where there may be no signal received for long segments of the 
propagation path* This behavior is expected, since randomness of 
atmospheric gradients and changing terrain features are common. The 
potential of several inversion layers existing is always there when 
the propagation path is great. 

In addition, for distances greater than 5 km, the received 
signal is fairly constant in level and the sound pressure does not 
follow the classical spherical divergence law. This variation from 
spherical spreading may be produced by the large number of multiple 
ray paths that are possible, with multiple ray arrival producing a 
mixture of phase that tends to produce a fixed sound pressure level. 

Since every sound propagation study in the long range is 
unique, the model was used to calculate the effect of changing a 
single parameter on the received signal. For instance, the source 
receiver geometry and the atmospheric refraction conditions were 
varied by selecting user parameters available from the program. The 
results of excess attenuation calculations were then compared for 
different standard days/nights. 

Figure 1 represents the total sound pressure level for the 
isotherraal-no wind condition for short detection distances of 200 
meters. This case illustrates the removal of refraction as an 
attenuation effect since the rays will travel in a straight line, 
with time of travel between equally spaced distances remaining the 
same. For low frequencies, especially 20 and 80 Hz, atmospheric 
absorption can be ignored and the curves illustrate the effect of 
spherical spreading and ground effects. 

The effects due to spherical spreading and atmospheric 
absorption were removed so that losses due to refraction and ground 
impedance could be examined more closely. Figure 2 examines the 
effect of isothermal atmospheres, where the excess attenuation is 
due to ground effects. Figure 2 shows that the model calculates the 
ground effect as a linear function of distance. 

Both atmospheric and wind refractive effects were investigated 
for the mid-latitude summer neutral profile for both the downwind 
and upwind cases, as seen in Figures 3 and 4. The excess attenuation 
is capped at 1 km and remains fixed for the entire range beyond 1 
km. For the upwind case, the cap starts at 2 km and the values 
remain fixed throughout the remaining ranges. One point should be 
made at this time; the values of excess attenuation for both cases 
are too low and refractive effects appear to be missing from 2 km 
onwards • 

The change in the meteorological profile to mid-latitude 
summer night is shown in Figures 3 and 6 for both wind directions. 
Again, the values are capped and the excess attenuation due to 
refraction is too low in value. 
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Consequently, there is a maximum distance beyond which the 
model should not be used. This distance is normally 1 km but can be 
extended to 2 km for atmospheric conditions that are unusually 
uniform. After 2 km, a model that uses instantaneous atmospheric 
readings to determine the velocity of sound profile should be used 
to calculate propagation losses. This latter model should use 
statistics determined by the defined topography and atmosphere to 
discuss variations in the received signal amplitude. 


CONCLUSION 

ADRPM*VII solves the detection problem even though detailed 
knowledge of temperature, humidity, variation in terrain features 
and wind gradients are not available to the user. Given these 
conditions, the model can give misleading information when compared 
to a model that performs ray tracing refraction based on accumulated 
local meteorological information. 

Perhaps a two model approach is required to solve the long 
range detection problem. ADRPM can be used for ranges below two 
kilometers where general meteorological conditions are approximated 
by readings at no more than two locations and terrain features are 
determined visually. Beyond two kilometers, a more elaborate model 
that is based on detailed atmospheric information would take over 
and perform the analysis. 
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Abstract 

One of the most important uses of acoustic propagation models lies in the area of 
detection and tracking of vehicles. Propagation models are used to compute transmission 
losses in performance prediction models and to analyze the results of past experiments. 
Vehicles can also provide the means for cost effective experiments to measure acoustic 
propagation conditions over significant ranges. In order to properly correlate the information 
provided by the experimental data and the propagation models, the following issues must be 
taken into consideration: 

• The phenomenology of the vehicle noise sources must be understood and 
characterized. 

• The vehicle's location or "ground truth" must be accurately reproduced and 
synchronized with the acoustic data. 

• Sufficient meteorological data must be collected to support the requirements of the 
propagation models. 

This paper treats the experimental procedures and instrumentation needed to carry 
out propagation experiments. Illustrative results are presented for two cases. First, a 
helicopter was used to measure propagation losses at a range of 1 to 10 Km. Second, a 
heavy diesel-powered vehicle was used to measure propagation losses in the 300 to 2200 m 
range. 

1. Introduction 

The development of acoustic propagation models has made significant advances in 
recent years resulting in accurate and practical propagation models such as those based on 
the Fast Field Program and the Parabolic Equation. Given sufficient meteorological data with 
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which to derive an accurate sound velocity profile, these programs model acoustic 
propagation losses quite accurately. Progress is also being made in the more difficult 
problem of modeling the effect of atmospheric turbulence on sound propagation. 

A very. significant sector of the uses for acoustic propagation models is in detection 
and tracking problems. In these systems, signals gathered by a microphone array are used to 
determine the location and track of a vehicle. Both air and ground vehicles are of importance 
in these applications. Acoustic propagation models play a very important role, being used to 
either predict performance under new conditions or to analyze the results of an experiment. 

In this paper we describe the methodology for the analysis of data involving vehicular 
sources and describe results obtained from two different tests: one, a long range experiment 
using a helicopter; the second, a mid range experiment, using a heavy, Diesel powered 
vehicle. 

2. Approach 

The essential elements necessary for the analysis of propagation data generated by 
vehicles are: a) a thorough understanding of the phenomenology of the vehicular sources, b) 
accurate positional data of the target vehicle's trajectory (ground truth data) and c) sufficient 
meteorological data to reconstruct the propagation conditions. 

2.1 Source Phenomenology 

In a test where the target vehicle is operating freely it is impractical to monitor the 
source strength continuously, therefore our knowledge of the source strength must be based 
on prior knowledge of the source's characteristics and whatever can be inferred by 
monitoring the observable parameters such as aspect angle or engine RPM. We will consider 
two types of vehicles, helicopters and heavy diesel powered vehicles. 

Helicopters provide an almost ideal source for long range propagation measurements. 
The noise generated by the main and tail rotors is loud and periodic with a relatively low 
fundamental frequency. In the spectral domain, helicopter signatures are characterized as 
families of narrow-band spectral lines. Fundamental frequencies of 10 to 30 Hz are typical. 
Source levels can reach 144 dB (re 20 micro-Pa in one Hz bands). Helicopters also operate 
at nearly constant blade rotational speed, as can be appreciated in a spectrogram (Figure 1), 
where the only frequency variations are those caused by the Doppler effect as the trajectory 
geometry changes. Strong aspect angle dependencies exist, both in the horizontal and the 
vertical planes (Figure 2). The amplitude of the rotor noise will also show a velocity 
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dependency proportional to the 12th power of the blade-tip Mach number (Figure 2). i-or a 
full treatment of helicopter noise characteristics, see Reference 1. 

Heavy Diesel-powered vehicles are easily detected at short to medium distances. Like 
helicopters, the spectral characteristics of vehicle noise are dominated by families of narrow- 
band harmonic components. Unlike helicopters, the frequency history of these components is 
highly variable. Rapid changes in engine RPM occur in response to operator actions, road 
conditions and gear changes. The amplitude of these narrow-band components is strongly 
dependent on engine load and RPM, as shown in Figure 3. From the sensor location, we 
must be content with observing only engine RPM. A good treatment of ground vehicle noise 
can be found in Reference 2. 

2.2 Vehicle Location Data 

Vehicle position data must be collected and synchronized with the acoustic data in 
order to measure propagation losses. In long range experiments or when the target is moving 
very fast, acoustic propagation delays must be accounted for. 

Helicopters and other aircraft can be tracked accurately with a radar system, if 
available. A more cost effective approach is to obtain tracking data from an Air Traffic Control 
facility, if the target is equipped with an ATC Beacon transponder. Such data can be obtained 
by prior arrangement with the local FAA facility. 

Ground targets can be tracked with an RF multilateration system, such as the Motorola 
Falcon Position Location System (PLS). This system is particularly convenient, since it allows 
tracking of multiple targets at a one Hz rate with digital data output. As an inexpensive 
alternative, the position of a ground target can be tracked by maintaining radio contact with 
one of the vehicle operators, who calls in "marks" as they go by pre-surveyed positions. 

Accurate ground truth is a necessity in these kinds of propagation experiments, but it 
need not be an inordinate expense if the proper procedures are worked out. 

2.3 Meteorological Data 

Meteorological data is a critical element of the propagation measurement, since it 
gives us the data necessary to understand the results of our experiment. 

The necessary meteorological information consists of sound velocity profiles, pressure 
and humidity. Sound velocity profiles are the most difficult to obtain. Traditional methods use 
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balloon sounding; SODAR devices are also being used at a limited number of sites. Tower 
measurements provide adequate data for short range experiments, and can fill the low 
altitude gap in the data provided by most balloon soundings. 

With progress being made in the modeling on the effects of turbulence on sound 
propagation, there is a need for more 'fine grained' measurements of the sound velocity 
profile. These gaps will have to be filled by more dense and frequent measurements of the 
lower atmosphere. 

3. Experiment Descriptions 

We will discuss the results of two propagation experiments using vehicular data. In the 
first a helicopter was tracked from a distance of 10 Km, in the second a diesel powered 
vehicle was tracked to a distance of 2.2 Km. 

3.1 Helicopter Test 

A test using a helicopter was made following a nearly radial trajectory starting at a 
distance of 1 0 Km. The helicopter flew at a speed of 1 85 Km and a height of 1 52 m. At the 
point of closest approach, it came within a distance of 500 m from the sensor site. The 
signature recorded by the sensors was shown in the form of a spectrogram in Figure 1 . The 
fourth harmonic was tracked automatically to extract frequency and amplitude data ( Figures 
4 ). Positional data was obtained with a radar tracking system and time-synchronized with the 
acoustic data. The constant speed trajectory allowed us to easily compensate for the 
propagation delays. 

Meteorological data consisted of a balloon sounding made approximately one hour 
before the test (Figure 5). The Fast Field Program was used to model propagation losses as a 
function of range, using the sound velocity profiles as input. A two parameter model of the 
surface acoustic impedance was used, with 300 Rayls of surface flow resistivity and porosity 
of 0.25. 

The results of the measured and modeled transmission losses (TL) are compared in 
Figure 6. Beyond a range of 4000 m the agreement between experimental and modeled data 
is quite good. The mean values of the TLs were very close. More important perhaps, the 
statistics of the variations with respect to their mean levels were also very close. It should be 
recognized that propagation losses will never be modeled beyond a certain level of precision 
and that a statistical description of propagation losses is the most realistic outcome given a 
limited amount of meteorological data. The statistics of signal and noise levels become 
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specially important in detection problems in order to predict the performance level of specific 
detection schemes. A study of the statistical properties of long range propagation losses 
appears to be a very promising area of research. 

The measured TL at ranges shorter than 4000m are higher than those predicted by the 
FFP. Some of the difference can be due to the directivity of the helicopter noise source, which 
reduces the effective source level as the elevation angle increases, however this effect is 
smaller than the observed discrepancy. At this point, we must attribute the differences to the 
inaccuracy of the sound velocity profile used in the FFP relative to the actual conditions at the 
time of the test. This result just reinforces the importance for accurate and timely 
meteorological data. 

3.2 Ground Vehicle Experiment. 

A short to mid range experiment was made using a heavy diesel powered vehicle. The 
vehicle operated on a road with a nearly radial trajectory starting at a range of 300m and 
finishing at a range of 2200 m. 

The vehicle signature as measured at the sensor location is shown in the form of a 
spectrogram in Figure 7. An automatic tracking program was used to extract the amplitude 
and frequency data corresponding to the 6th engine harmonic or Engine Firing Rate; this 
information is shown in Figure 8. 

Lack of sound velocity profile data forced us to model the SVP as that of a 'neutral' 
atmosphere, that is, a profile matching the nominal atmospheric lapse rate. The neutral 
atmosphere profile was used as an input to the FFP, producing the TL curve shown in Figure 
9, along with the measured TL. The match between the measured and modeled TLs is good 
at short ranges, but they start to diverge at longer ranges. However, a simplistic model which 
assumes spherical spreading plus a fairly high absorption term ( 0.0045 dB / meter) 
produced an excellent fit to the measured data. We hypothesized that the difference could, in 
part, be explained by variations in the engine RPM and/or engine load. The noise of heavy 
diesel powered vehicles is directly affected by engine load and RPM. An attempt was made 
to compensate for the effect of RPM. This is an imperfect approach, since we should 
compensate for both the RPM and load, however we do not know of any practical way of 
inferring load at long distances. The incremental sound pressure level relative to the best 
fitting model was plotted against the incremental frequency relative to 80 Hz. This result is 
shown in Figure 10, and shows a clear dependency between SPL and frequency. The SPL 
figures were then adjusted to a constant 80 Hz (SPL was adjusted downward if the frequency 
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was more than 80 Hz, upwards if it was less than 80 Hz), producing the curve shown in 
Figure 1 1 . The corrected TL curve shows a better agreement with the computational models. 
Some of the extreme variations in TL have also been reduced as a result of the 
compensation procedure. 

4. Conclusions 

Two experiments involving a ground vehicle and an aircraft have been analyzed with 
the help of the Fast Field Program, one of the state of the art acoustic propagation models. 

By making use of meteorological data as an input to the Fast Field Program and knowledge 
about the source phenomenology of the vehicles, we were able to obtain a good match 
between the measured and predicted transmission losses. These results are encouraging 
and underscore the importance of thoroughly characterizing vehicular sources and of 
obtaining fine grained meteorological data. 

The development of computational models of sound propagation have made dramatic 
advances in recent years, and their need becomes the driving requirement for data collection 
in many field experiments. 
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Figure 1 Spectrogram of helicopter signature, recorded during a nearly 

inbound-radial run. 
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Figure 2 Helicopter noise sources are highly directional and strongly 
dependent on blade tip velocity (and therefore on helicopter forward speed as 

well). 
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Figure 3 The noise generated by a Diesel-powered ground vehicle is strongly 
dependent on engine speed (RPM) and load. 
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Figure 4 Amplitude and frequency of the fourth harmonic of the helicopter 

noise signature as a function of time. 







Figure 5 Profiles of temperature and wind velocity prior to the propagation 

experiment. 
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Long Range Experiment: TL vs Range 
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Figure 6 Measured and modeled (using the FFP) transmission losses. 
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FREQUENCY 


Figure 7 Spectrogram of the acoustic signature of a heavy Diesel-powered 

vehicle on level ground. 
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Figure 8 Amplitude and frequency of the main component of the engine noise 

signature as a function of time. 
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Figure 9 Measured and modeled (using the FFP) transmission losses during 

outward-bound run. 
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Figure 10 Measured Sound Pressure Level shows a strong dependency on 
engine RPM. Engine load was not directly observable in this experiment. 
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Figure 11 Measured transmission loss, after correcting for the effects of 

engine RPM variations. 
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SCATTERING MEASUREMENTS ON NATURAL AND MODEL TREES ‘ ' 


James C. Rogers and Sung M. Lee 
Michigan Technological University 


SUMMARY 


The acoustical back scattering from a simple scale model of a 
tree has been experimentally measured. The model consisted of a 
trunk and six limbs, each with 4 branches; no foliage or twigs were 
included. The data from the anechoic chamber measurements were 
then mathematically combined to construct the effective back 
scattering from groups of trees. Also, initial measurements have 
been conducted out-of-doors on a single tree in an open field 
in order to characterize its acoustic scattering as a function 
of azimuth angle. These measurements were performed in the spring, 
prior to leaf development. The data support a statistical 
model of forest scattering; the scattered signal spectrum is 
highly irregular but with a remarkable general resemblance to 
the incident signal spectrum. Also, the scattered signal's spectra 
showed little dependence upon scattering angle. 


INTRODUCTION 


Acoustic scattering in forests has often been studied in the 
context of sound which propagates through forests and thereby suffers 
attenuation. This attenuation is attractive to those who might 
consider the acoustic screening effects of forested areas. Thus, 
sound propagation in forested areas has been considered by many 
researchers (ref. 1, 2, 3). At least five factors contribute to the 
attenuation of sound propagating in forests: spherical spreading, 
atmospheric absorption, foliage absorption, ground loss, and 
scattering. It appears that scattering is a significant factor in 
sound attenuation at the middle frequency range (ref. 4, 5) . The 
approach toward studying scattering that we use here is to focus on 
scattering alone and to particularly include back scattering. In 
this way only the scattered signal is measured whereas in 
traditional measurements of attenuation through forests both 
scattered signals and direct signals are present. In this case it 
is quite difficult to separate the scattered component from the 
considerably stronger direct signal component. Since forests are 
made up of many single trees, back scattering from forests can be 
considered using single tree scattering processes and extending this 
to the aggregate effects of many trees. 
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MODEL TREE BACK SCATTERING MEASUREMENTS IN AN ANECHOIC CHAMBER 


Purpose of Measurements 

There are several advantages to making back scattering 
measurements on a simple model tree in an anechoic chamber. If an 
asymmetrical tree is used the scattered signal will be different for 
different azimuth angles and synthetic "forest scattering" data can 
be generated by using an ensemble of these model scattering 
measurements. Also, it is possible to observe scattering as a 
function of increasing scattering angle: zero degrees for back 
scattering and one hundred eighty degrees for forward scattering. 

The results reported here do not use this capability however. If 
the tree is elevated on a pedestal the effects of ground reflections 
are removed, something that is not possible with a natural tree. 
Finally, the measurements are quite repeatable with no effects from 
meteorological influences. 


Description of Model "Tree" 

A tree silhouette was selected that approximately simulates 
that of a tree in a northern hardwood forest. The basic structure 
is a trunk, several limbs and a large number of branches as 
described earlier in Rogers et. al. (ref. 6). In order to utilize a 
simple construction technique and to facilitate theoretical analysis 
(something not yet completed) , the cylinder shape was used as a 
basic structure element in our model. Hard wooden cylinders of 
three diameters were used for fabricating the three basic elements: 
the trunk, the limbs, and the branches. A single wooden cylinder 
that is several wavelengths long provides an effective back 
scattering element with a structured scattering pattern (ref. 7-) . 
Figure 1 shows a sketch (not to scale) of the tree and lists the 
dimensions and numbers of the components. The limbs were randomly 
distributed around the perimeter of the trunk and were spaced at 
irregular intervals along its length. The branches were similarly 
placed on the limbs. The effects of leaves and twigs were ignored. 
We believe that these will not give significant back scattering 
contributions in the low and mid frequency range studied. 


Back Scattering Measurements and Results 

The basic arrangement of the speaker source and receiving 
microphone in relation to the tree are shown in Figure 1. A single 
pulse was applied to the speaker through an amplifier, received by 
the microphone as a "direct" wave, and again received by the 
microphone as a back scattered signal. 

Our small speaker (with a hemispherical cone approximately 
0.02m in diameter) did not radiate a great deal of energy and 
several techniques were used to ensure an adequate signal to noise 
ratio. The tree was removed from its stand and a "constant 
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background" measurement was made by coherently time averaging 
several pulse events; this gave a reliable estimate of the signal 
which regularly existed in the portion of a time record occupied by 
the desired scattered signal. This signal was subtracted from all 
scattering records. Also, coherent time averaging was used in all 
scattering measurements to reduce the effects of random noise. 

Using the known geometry it is possible to construct a time window 
in which the scattered events will appear; such a window was used to 
exclude all signal outside of the desired scattering events. 

Finally, the useful spectral content of the source was judged to be 
from approximately 1 kHz to over 10 kHz. A filter was applied to 
the scattered signal to permit only those frequencies in the 
analyzed data. 

Figure 2 compares the signal back scattered from the trunk 
alone, after it has been processed as described above and amplified 
by a factor of approximately 30, with the "direct" signal. There is 
a high degree of similarity between the signals as would be expected 
for back scattering from a single cylinder. The scattered signal is 
considerably more complex after the limbs and branches are added to 
the trunk. 

After the tree was assembled, twenty four separate back 
scattering measurements were made. For each the tree was rotated 15 
degrees about its vertical axis. The non symmetrical nature of the 
silhouette produced 24 unique scattering records which were then 
treated as the back scattering from 24 separate trees. 


SYNTHESIS OF BACK SCATTERING FROM A GROUP OF MODEL TREES 


Eighteen of the unique back scattering records were used to 
synthesize the scattering one would measure from a grove of 18 
trees. Figure 3 shows a plan view of the grove, to scale, where the 
distance between the source, the microphone, and the first tree in 
the grove is indicated. Each original time domain measured back 
scattered signal was amplitude scaled by a (1/distance) factor to 
account for the round trip distance from the source to the tree and 
back to the microphone. Also, each original signal was time delayed 
by an amount proportional to the round trip distance. Finally all 
eighteen time domain records were added to simulate the signal back 
scattered from a grove of trees. Figure 4a shows the composite 
time domain signal. Two features are apparent; several individual 
tree scattering events are seen and as time increases in the figure 
the signal amplitude diminishes in accordance with the (1/distance) 
spreading factor. Figure 4b presents the spectrum of the composite 
scattered signal and compares it with the spectrum of a single 
direct pulse. 

Over the useful bandwidth of the signal shown in the figure, 
approximately 1 kHz to 14 kHz, there is a close resemblance between 
the average spectrum of the scattered signal and that of the direct 
signal. There is approximately a 30 dB level difference between the 
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two spectra and a very irregular character to the scattered signal 
spectrum as would be expected for a random combination of similar 
signals. This random spectral appearance is observed even though 
the model tree did not have a broad distribution of sizes in its 
structure (in fact, only three different cylinder diameters and 
lengths were used) . We conclude, by comparing the spectrum of the 
scattered signal from the synthetic grove with that which is 
produced by an actual forest (not shown here) that a high degree of 
realism has been achieved with a relatively small number of "trees" 
in the grove. One further comment about the synthetic scattering 
record should be made: since we combined individual records of 

sound scattered from individual trees, we have not allowed multiple 
scattering between trees. However, each individual tree record 
naturally incorporates multiple scattering among elements of the 
tree such as the trunk, limbs and branches. This scattering is 
probably considerably more important than that between individual 
trees . 


SCATTERING MEASUREMENTS ON A SINGLE TREE IN A FIELD 


Purpose and Measurement Arrangement 

Scattering within a forest is a complex process; the presence 
of a large number of individual scattering trees with a wide spatial 
distribution precludes the study of the process at the level of the 
individual tree. We have therefore selected an isolated tree 
located in a uniformly flat grassy field for a series of scattering 
measurements. Both back scattering and scattering at angles up to 
165 degrees from back scattering have been measured. This 
arrangement permits use of an impulsive source which is desired for 
separating the scattered signal from the signal which travels 
directly from the source to the microphone. The source was a simple 
mechanical device with a barrel and firing pin. It was machined to 
accept shot shell primers, Winchester part # 209, which are 
detonated by striking the firing pin. 

Figure 5 is a plan view of the measurement arrangement. The 
source was located at a fixed point 30 meters from the center of the 
tree and the measurement microphone was located a distance of 15 
meters from the tree at a series of points separated by 15 degrees 
of azimuth. A reference microphone was situated along a line 
between the source and the tree and 5 meters from the source. The 
source and measurement microphones were at fixed heights of 1.15 and 
1.10 meters respectively. Bruel and Kjaer microphones, type 4155, 
were used on type 2330 sound level meters for both the reference and 
measurement microphones. Typical peak direct wave sound levels 
measured by the reference and measurement sound level meters were 
137 dB and 120 dB respectively. At each location three separate 
shots were fired and the data recorded on a 4 channel digital audio 
tape recorder with a uniform frequency response from 0 to 10 kHz and 
a dynamic range of 84 dB. 
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Ground Impedance Measurements 

It is not possible to directly remove the effects of ground 
reflections from the measurements since the source-ground-tree 
geometry is quite variable over all of the tree components such as 
the trunk, limbs, and branches. A level difference measurement was 
made between two microphones; one was placed at the ground level 
and the second was elevated 1.15 meters directly over the first. 

The shot source was located 5 meters away from the pair at an 
elevation of 1.15 meters. The ratio of the elevated microphone 
power spectrum to the ground level microphone power spectrum 
produced a differential spectrum or "transfer function" magnitude 
characteristic of the interference process between the direct and 
the reflected wave as described in ref. 4. Using the experimental 
data found in Donato (ref. 8) and the fact that the real and 
imaginary parts of the ground impedance are observed to vary 
approximately as the inverse square root of the frequency, a good 
fit was found for our experimentally determined differential 
spectrum with a theoretically predicted differential spectrum. The 
fit was better at frequencies below 1000 Hz but quite acceptable 
above that frequency too. We thus have a good estimate of the 
ground impedance for the field surrounding the tree. A 
representative value for the magnitude of the ground reflection 
coefficient at 300 Hz is about 0.8. 

The initial scattering investigation sought to minimize the 
variability of all effects except the azimuth angle which was varied 
in 15 degree increments as shown earlier. Thus, although the 
precise effect of the ground reflections on the " insonif ication 
function" for the tree is not known, the source-tree geometry was 
fixed for all of the measurements. Also, the measurement microphone 
was always maintained at a fixed distance of 15 meters from the 
tree. The impulsive source proved to be quite repeatable but to 
reduce the effects of random noise and source variability somewhat, 
each scattering measurement reported here is the average of three 
power spectra from three separate measurements. The temperature was 
approximately 78 degrees F at 1 meter elevation and the wind varied 
in strength from 0.5 to about 1.5 meters per second. 


Results of Scattering Measurements 

Figure 6 shows the back scattered (0 degree azimuth angle) 
signal spectrum. The general shape is characteristic of the source 
alone in an anechoic environment (without any ground effect 
present) . The spectrum, which is the average of three power 
spectra, is highly irregular in the same manner as that previously 
observed in Figure 4b for the synthetic grove of trees. The maple 
tree used for the outdoor experiment had multiple trunks with 
dozens of limbs and branches. Since the measurement was made early 
in the spring there were no leaves on the tree. One can estimate 
the signal to noise level by examining the background noise spectrum 
(the average of three noise power spectra) which is also shown on 
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the figure. There is good signal to noise (about 15 dB to 25 dB) 
over an approximate frequency range from 0.5 kHz to 9.5 kHz. 

The scattered signal spectra at azimuth angles of 45, 90, and 
150 degrees are shown in Figures 7, 8, and 9 respectively. Also, 
the scattered signal at 0 degrees, from Figure 6, is shown in these 
figures. A principal feature of these plots is that the scattered 
signal spectrum level varies only a small amount with azimuth angle 
for the frequency range of the measurements. 


DISCUSSION AND CONCLUSIONS 


Close examination of Figures 7 , 8 , and 9 shows that the 
scattered signal is remarkably similar over all angles of 
scattering. The insonif ication of the tree was the same in all 
cases since the source-tree geometry remained constant so comparison 
between figures examines only the effects of varying the scattering 
angle. Only at 150 degrees (Figure 7) does there appear to be a 
significant variation from the scattered signal at 0 degrees. This 
deviation is seen in frequencies from approximately 700 Hz to 1200 
Hz. The rather narrow frequency range of this feature is 
perplexing. The scattering record for 165 degrees has been examined 
in this frequency range. It too shows reduced signal levels over 
approximately the same frequency range. However, the scattering 
record at 135 degrees does not display this feature. 

The measurements presented here are part of a continuing 
investigation of scattering by forests. Foliage effects are to be 
included and additional low frequency data are required. Also, 
scattering models for trees and forests are required. These should 
adequately treat azimuth angle, frequency effects and address the 
problem of ground effects. Finally, the influence of tree variety 
should be considered on forest scattering. 
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Amplitude 


AN ARTIFICIAL TREE 


TREE SPECIFICATIONS: 


1 - TRUNK. 

6 - LIMBS. 

24 - BRANCHES 


l.Om x 0.050m Ola, 
0.4m x 0.025m Ola. 
0.2m x 0.012m Dla, 


ASYMMETRICAL ARRANGEMENT 
GEOMETRY: 

Source - Receiver - 1.05 m 
Receiver - Tree - 0.87 m 


Sketch of model tree (not to scale) with element 
dimensions and physical arrangement for back 
scattering measurements in an anechoic chamber. 



TIME, MILLISECONDS 

The direct signal (smooth) and the back scattered 
signal from the trunk alone (noisy and irregular) 
after amplification by a factor of 30. 







Figure 3. Arrangement of individual trees in synthetic grove 

of trees (to scale) . The source and microphone 
positions are also to scale. Each tree position 
contributes an individual time domain scattering 
event to the synthetic back scattering record as 
described in the test. 
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AVERAGE OF THREE SPECTRA AT 45 DEGREES 



e 7. Single tree scattering at 45 degrees, solid line, 

compared with back scattered spectrum, dashed 
line. 


AVERAGE OF THREE SPECTRA AT 90 DEGREES 



FREQUENCY, HZ 

Single tree scattering at 90 degrees, solid line, 
compared with the back scattered spectrum, dashed 


Figure 8. 
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Abstract 

Using a Fast Fourier integration method and a global matrix method for solution of the boundary 
condition equations at all interfaces simultaneously, a useful tool for predicting acoustic propagation 
in a stratified fluid over a stratified porous-elastic solid has been developed. The model for the solid is 
a modified Biot-Stoll model incorporating four parameters describing the pore structure corresponding 
to the Rayleigh- Attenborough rigid-porous structure model. 

The method is also compared to another Fast Fourier code (CERL-FFP) which models the ground 
as an impedance surface under a horizontally stratified air. Agreement with the CERL FFP is good. 

The effects on sound propagation of a combination of ground elasticity, complex ground structure, 
and atmospheric conditions are demonstrated by theoretical results over a snow layer, and experimental 
results over a model ground surface. 


Introduction 

The ground has conventionally been modelled for outdoor sound propagation as either an impedance 
surface or a rigid-porous structure. These approaches have both been highly productive in the case of 
high density materials. However in reality the ground is poro-elastic. Ground surfaces have hitherto been 
modelled as such when the interest has been in acoustic to seismic coupling, but there has been little 
interest in porous-elastic ground models in propagation in the air. For some outdoor ground surfaces 
(such as snow or forest floors for example) the bulk density of the material is low enough for seismic 
effects to become important for sound propagation over the surface at some frequencies. 

In this paper an FFP propagation model is used to calculate sound pressure levels over a porous-elastic 
ground surface. The model’s predictions are compared to the predictions of other propagation models for 
the high density, high seismic velocity rigid-porous limit of the porous-elastic ground model. The effects 
on acoustic propagation of the elasticity of various ground surfaces is then shown by comparison to the 
rigid frame limit. Using a multiply layered fluid atmosphere the combined effects of meteorology and 
ground elasticity are examined. 

The Biot-Stoll poro-elastic model 

The ground model used in this investigation was a modified Biot-Stoll Poro- Elastic model[l,2,3]. Propa- 
gation within the material is via three different modes; a fast wave, equivalent to the seismic P wave: a 
slow wave equivalent to the pore wave in the Rayleigh Attenborough rigid-porous model[4]: and a shear 
wave equivalent to the seismic S wave. Each wavetype causes vibration in both the solid material and 
the pore fluid. Attenuation of all three wavetypes is predicted by the theory due to viscous losses on the 
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pore walls, though it has been shown to underpredict the attenuation in real materials because other loss 
mechanisms are not taken into account. Hence an extra attenuation is added as an imaginary part of the 
fast and shear propagation constants. 


The Fast Fourier Method 

If one applies a Hankel transform in range to the Helmholtz equation one obtains the depth separated 
wave equation: 


where, for a point source, 


(j[^ + (k 2 - /&(*)) ^ r (k, z) = h(k, z), 


h(k,z) = —6(2 - z Q ) 


Cl) 


( 2 ) 


Solutions to this equation are depth dependent only and are equivalent to solutions to the wave equation 
for continuous plane wave incidence. In order to obtain a range dependent solution one must obtain 
depth dependent solutions to the depth separated wave equation, and then perform the inverse Hankel 
Transform on the solution to equation 1 , which is in terms of horizontal wavenumber. 

The exact range dependent solution is in the form; 


r 00 

J kh=0 


!k h 

where T is the depth dependent Greens function. 

A large argument approximation to the Bessel function [5] is: 
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This approximation together with the replacement of the integration by a finite sum gives the approximate 
equation for ^ : 
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This Fourier series approximation can then be improved by corrections to allow for the truncation of 
the integral to infinity to a finite wavenumber, ^( maI ), and the avoidance of pole(s) on the real axis[6], 
which together lead to inaccuracies and oscillations in the result, to give 
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and 5’ is an approximation to the sum, 
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The Environment 




The environment is assumed to be range independent and to consist of a fluid (air) upper half-space 
overlying a set of horizontal fluid (air) layers of differing sound speeds and densities. The lowest of these 
fluid layers is in contact with a ground made up of a set of horizontal elastic porous layers under which 
is an elastic porous half-space. The number of layers in either fluid or ground can be set to zero. 


The Depth Dependent Green’s function 


The depth dependent Green’s function T must be solved for the above environment. 

In a fluid layer containing a spherical point source the depth dependent Green’s function is 


T = 


( 10 ) 


The iZj are calculated by solution of the boundary condition equations at the interfaces. 

In the porous elastic medium there are three scalar displacement potentials describing propagation 
in the fluid, 
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Jo 
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3*1 is the longitudinal displacement potential in the solid, $2 is the longitudinal displacement potential 
in the pore fluid, $ 3 is the transverse displacement potential in the solid, to which the fluid transverse 
displacement potential is directly proportional. 

In a porous-elastic layer, bounded by interfaces at depths d\ and d 2) in the absence of a source, the 
$iS at a depth z are given by. 


*1 = A n e< z ~ d ^ + + A 2l e^ dl ^ + A 2 ^~ z ^ , (14) 

*2 = m! (A n e^ z - d ^ + + m 2 (A 2i e< z ~ d ^ + A 2T e‘ (, * a ~* )ft ) , (15) 

$ 3 = A 3l e^ z ~ d ^ 3 + A 3 ^~ z ^ 3 , (16) 

The m.i are the ratios of the amplitude in the solid and pore fluid for each wavetype, and the /3, = 
( k 2 - fc 2 ) 1 / 2 , where the ki are propagation constants, and k is the horizontal wavenumber. The depth 
dependent Green’s function T for a desired output parameter in the fluid is a function of the The Aj 
are calculated by solution of the boundary condition equations. 

Boundary conditions 

Boundary condition equations in cylindrical polar coordinates (r,0,z) are needed. However the axisym- 
metric nature of the problem considered here means that there is no 6 dependence. 

At boundaries between two fluid layers the two boundary conditions are 

1 . continuity of pressure, 

2. continuity of normal particle displacement, 

At the interface between the fluid and the porous elastic medium there are four boundary conditions, 
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1. continuity of total normal stress [3], 

2. continuity of normal displacement [3], 

3. continuity of fluid pressure[7], 

4. continuity of tangential stress, 

Six boundary conditions are required at each interface between porous-elastic layers. These boundary 
conditions are the four above and two others, 

5. continuity of normal relative fluid displacement, and 

6. continuity of tangential frame displacement. 

The range dependent parts of the boundary condition equations are identical on each side of the boundary, 
therefore only the depth dependent Green’s functions of the boundary conditions need to be equated. 

The boundary condition equations are solved simultanaeously for the A j and i?j at all interfaces. The 
depth dependent Green’s function is then calculated for the desired output parameter (sound pressure 
level, frame displacement, etc). The range dependent solution is calculated using the FFP method 
described above. 

Comparison to other propagation models 

For propagation above a rigid-porous halfspace the model compares well with other propagation mod- 
els, such as the CERL-FFP (see figure 1)[8], and Attenborough, Hayek, and Lawther’s ’exact’ analytic 
model(see figure 2)[7] . Above an extended reaction rigid-porous layer over a non-porous backing agree- 
ment with Nicholas-Berry and Daigle’s propagation model is good for a wide variety of model surfaces(see 
figures 3 and 4) [9]. Source and receiver heights are 0.5 and 0.3 metres respectively. 

Effects of ground surface elasticity on sound propagation 

The largest effects of ground elasticity on sound propagation over it are likely to be where the bulk density 
of the ground surface is small. The most common ground cover where this is so is a snow layer. Measured 
normal surface impedance over snow cover sometimes shows low frequency peaks [10,11], These could be 
interpreted as seismic resonances in a snow layer. Figure 5 shows the predicted excess attenuation over an 
8cm thick snow layer overlying a rigid nonporous halfspace at twenty metres range, using a rigid-porous 
model, and porous-elastic model. The pore structure and elastic parameters are taken calculated from 
Sommerfeld[12], Johnson[ll] Ishida[l0] and Attenborough and Buser[l3]. A resonant effect can clearly 
be seen at about 810Hz in the porous elastic model output which is not present for the rigid-porous 
model. Figure 6 shows the predicted excess attenuation over the same snow layer at 810Hz as a function 
of range. This figure demontrates that at this frequency a seismic resonance in the snow layer leads to 
an apparent hardening of the snow surface at a short range, leading to less attenuation due to ground 
absorption. The behaviour at longer ranges shows that away from the source the attenuation due to 
the ground is unaffected by the elastic effects, but the signal amplitude is increased due to the reduced 
attenuation near to the source. 

Combined effects of elasticity and atmospheric sound velocity gradients 

Continuous sound velocity gradients can be modelled by thin homogeneous layers as long as the layer 
thickness is much less than the wavelength of the sound [14], In figure 7 the combined effect of the 
logarithmic downward refracting sound velocity gradient(roughness length 5.10 _3 metres, temperature 
difference between ground and 4.0 metres 1° Centigrade) and an elastic surface are shown. The difference 
between elastic and rigid models remains approximately the same as for no gradient. 
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Comparison with experiments 

In order to test the validity of this porous-elastic propagation model, measurements of the level 
difference between two vertically separated microphones were made over a thin (4cm) layer of low density 
foam material. The foam was attached to a non-porous concrete surface. A point noise source was 
suspended over the foam surface. 

The elastic and porous parameters of the foam were separately measured using non-acoustic tech- 
niques. The measured level difference was compared to the level difference predicted using both rigid 
and poro-elastic models. The results are shown in figure 8. The geometry used for this figure was source 
height 0.2 metres, receiver heights 0.01 and 0.2 metres, and range 0.4 metres. The results show a better 
agreement with the elastic model than with the rigid model. 

Conclusions 

An FFP model for propagation over porous-elastic surfaces has been developed. It has been shown that 
in the rigid frame limit it agrees well with other propagation models. For sound propagation over low 
bulk density layered materials it has been shown that ground elasticity can have a substantial effect on 
received sound pressure levels for both real and theoretical results. 
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Table 1: Material parameters used in the prediction of excess attenuation 


Parameter 

Unit 

Rigid-porous 

Halfspace 

Rigid-porous 

layer 

Snow 

layer 

Foam 

layer 

Flow resistivity a 

MKSraylsm” 1 

100000 

10000 

15900 

18400 

Porosity Q 

- 

0.3 

0.3 

0.804 

0.97 

Pore shape factor ratio s p 

- 

0.5 

0.5 

0.5 

6.5 

Grain shape factor n* 

- 

0.5 

0.5 

0.5 

63.8 

Bulk density 

kgmr 3 

- 

* 

184.0 

32.0 

P-wave velocity v p 

771S” 1 

- 

- 

130.0 

79.0 

S-wave velocity v 3 

ms” 1 

- 

- 

90.0 

56.0 

Q(v)/»( v) 

- 

- 

- 

0.05 

0.085 

Grain bulk modulus K r 

Nm~ 2 

- 

- 

1.10® 

1.10 10 

Layer depth 

m 

- 

0.1 

0.08 

0.04 
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Ranqe(m) 

Figure 1 Comparison of CERL FFP and FFLAGS for propagation over a 
rigid-porous halfspace in a 12 layered atmosphere. Frequency=50Hz. 



Figure 2 Comparison of FFLAGS to the predictions of Attenborough, Hayek, 
and Lawther’s exact extended reaction model for propagation over an 
extended reaction rigid-porous halfspace. Range=20 metres. 
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Figure 3 Comparison of Nicholas Berry and Daigle’s model for predicted ex- 
cess attenuation over an extended reaction rigid-backed layer at 100Hz. 



100 500 1000 5000 

Frequency(Hz) 

Figure 4 Comparison of FFLAGS to Nicholas Berry and Daigle’s model for 
predicted excess attenuation over an extended reaction rigid backed 
layer at 20 metres. 
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Frequency(Hz) 
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Figure 5 Predicted excess attenuation over thin snow layer using FFLAGS 
with rigid-porous parameters, and porous-elastic parameters. Range 
20 metreB. 



Figure 6 Predicted excess attenuation over thin snow layer at 810Hz. Show 
ing difference between rigid and elastic model at this frequency. 






Frequency(Hz) 

Figure 7 Predicted excess attenuation over thin snow layer at 810Hz in the 
presence of a downward refracting logarithmic sound velocity gradient. 
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Frequency (Hz) 


Figure 8 Measured Level difference between two vertically separated micro- 
phones (heights 0.2 and 0.01 metres) at a range of 0.4 metres from a 
point source over a thin (4cm) rigidly backed foam surface. Compared 
to the predicted level difference using rigid-porous and porous-elastic 
models, calculated from direct measurement of the pore and elastic 
parameters of the foam. 
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N91-16693 


William L. Willshire, Jr. 
NASA Langley Research Center 
Hampton, VA 23665 

and 


Donald P. Garber 

Lockheed Engineering & Sciences Company 
Hampton, VA 23666 


Development of the advanced turboprop has led to concerns about en route 
noise. Advanced turboprops generate low-frequency, periodic noise signatures at 
relatively high levels. As demonstrated in a flight test of NASA LeRC's Propfan Test 
Assessment (PTA) airplane in Alabama in October 1987, the noise of an advanced 
turboprop operating at cruise altitudes can be audible on the ground. The assessment 
of the en route noise issue is difficult due to the variability in received noise levels 
caused by atmospheric propagation and the uncertainty in predicting community 
response to the relatively low-level en route noise, as compared to noise associated 
with airport operations. 

The En Route Noise Test was designed to address the atmospheric propagation 
of advanced turboprop noise from cruise altitudes and consisted of measuring the 
noise of an advance turboprop at cruise in close proximity to the turboprop and on 
the ground. Measured and predicted ground noise levels will be presented in this 
paper. Participants in the En Route Noise Test were NASA LeRC, the FAA, and 
NASA LaRC. 


EXPERIMENT DESCRIPTION 

The test airplane was NASA LeRC's PTA airplane which has a 2.7 m (9 ft) 
diameter, eight-bladed, tractor-configured advanced turboprop mounted on its left 
wing. The test airplane was instrumented to measure the near-field turboprop noise 
levels, as well as, engine and other pertinent parameters. During the microphone 
array flyovers, the test airplane was tracked with a C-band beacon. 

The En Route Noise Test was performed at the White Sands Missile Range in 
New Mexico in April 1989. Eighty-eight PTA airplane passes or runs over the ground 
microphone array were recorded. The array was an eight element, linear 
microphone array with an inter-element spacing of 122 m (400 ft). The completed 


127 



test matrix is illustrated in Table I. The majority of the runs were performed at 
altitudes of 4.6 and 9.2 km at a tangential tip speed of 240 meters per second (bpf of 
, 226 Hz) and a aioihinal power setting of 90 percent. Seventeen runs were flown at 

other tip speeds in the range of 190 to 260 m/s. Meteorological profiles were 
measured during the flyovers from ground level up to 12 km. 


PTA 

SPEED, M 

ALTITUDE, km AGL \ 

.6 

2.7 

4.6 

9.2 

.5 

4 

4 

23 


.7 



19 

32 

.77 




6 


Table I. Completed test matrix. 


DATA ANALYSIS 

The basic analysis used in the results presented in this paper is ensemble- 
average time histories^ . Data from the eight ground mounted digital microphones 
(SR=2344 sps) are high-pass filtered at 80 Hz and then shifted based on the airplane 
ground speed to give all eight individual microphone time histories a common 
source emission time base. Each individual microphone time history consists of a 
series of 1 /2-second root mean square pressure levels. The shifted time histories are 
then averaged. An ensemble-average time history has less variability than a single 
microphone time history and increased statistical confidence. 

RESULTS 

Data Variability .- To investigate long-term, between day, data variability, peak 
Overall Sound Pressure Level (OASPL) for each run was calculated from the 
ensemble-average time histories. No corrections were applied for deviations from a 
nominal flight path and no runs were rejected. The peak levels were averaged for 
like test conditions on a daily basis. Results are given in Table II for the 4.6 and 9.2 
km runs with a tip speed of 240 m/s. Average OASPL valves ranged from 61 to 
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75 dB. In general, the repeatability within a particular test day was good. The 
standard deviation of the average OASPL of the 11 similar rims which were flown in 
a 90 minute period during the 6th test session was .7 dB. However, the average 
levels for the same test condition varied from day to day. For the 9.2 km, .7 M test 
condition, there was an 11 dB difference in average levels across days. For the same 
test condition and runs, a boom microphone on the PTA aircraft exhibited a range of 
3 dB in the blade passage frequency noise level. The ground measured OASPLs were 
dominated by the blade passage frequency sound pressure level. This indicates, as 
expected, that the variability observed in the ground measurements is largely due to 
propagation. Another observation is that on the 3rd day the .5 M, 4.6 km average 
levels are greater by 3 dB than the .7 M, 4.6 km average levels. 


TEST SESSION 


TEST 

KEY 

1 

2 

3 

4 

5 

6 

7 

8 

CONDITION 

AVG, dB 

60.7 

69.0 

60.7 

65.1 



67.4 

72.1 

9.2 km, .7 M 

a, dB 

1.4 

.6 

.2 

.9 



2.2 

.8 


No. 

2 

4 

4 

4 



3 

4 


AVG, dB 

74.8 

72.6 

67.6 

69.7 

74.8 



74.0 

4.6 km, .7 M 

o, dB 

1.8 

.5 

1.3 

.9 

1.5 



1.9 


No. 

2 

2 

4 

4 

3 



4 


AVG, dB 

72.2 


70.6 

70.2 

74.7 

74.3 



4.6 km, .5 M 

o, dB 

.6 


.9 

.2 

.1 

.7 




No. 

2 


4 

3 

2 

11 




Table n. Averaged ensemble-average time history peak OASPL. 


Comparison To Ray Tracing .- Figure 1 is a comparison of a ray tracing predicted time 
history to an ensemble-average time history. The measured data are from a 9.2 km, 

.7 M run with a tip speed of 240 m/s. Included in the figure are the ensemble-average 
80 percent confidence bounds. The acoustic source used in the ray tracing 
propagation model was an ANOPP^ prediction based on nominally measured 
advanced turboprop operating conditions. An amplitude correction was applied to 
the predicted source levels for each run type based on the difference between a 
predicted and measured boom microphone amplitude for each run. A radiosonde 
weather profile was used in the two-dimensional ray tracing model which 
incorporates the effect of the wind by calculating an effective sound speed which 
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includes the component of the assumed horizontal wind in the vertical plane 
containing the airplane and the receiver. This comparison between measured and 
predicted time histories is fair. The peak level is overpredicted by 4 dB, and there is a 
small time shift in the prediction. 



Figure 1. Ray trace prediction compared to ensemble-average time history. 


Prediction Error .- To illustrate prediction error, the measured versus predicted 
overall SPLs for the 9.2 km, .7m; 4.6km, .7 m, and the 4.6 km, .5 m runs are plotted, 
respectively, in figures 2a through 2c. The dashed line in the plots is the perfect 
agreement line. The middle solid line is a regression line, and the 80% confidence 
levels about the regression line are represented by the two remaining solid lines. For 
the first two test conditions, 9.2 km and 4.6 km with a .7 m, the perfect agreement 
line falls within the 80% regression confidence bounds. There is an approximately 2 
dB underpredicted basis in the 4.6 km, .5 m results. The reason for the basis is not 
currently known. The procedure for estimating source levels is being carefully 
reviewed. In general, the agreement between measurement and prediction is judged 
to be good. 
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asured levels, 
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Figure 2b. Prediction error result for Mach .7, 15,000 ft. altitude condition. 









Figure 2c. Prediction error result for Mach .5, 15,000 ft. altitude condition. 

SUMMARY 

A flight experiment was conducted to investigate the propagation of advanced 
turboprop noise from cruise altitudes. The experiment was designed to use ensemble 
averaging and to measure weather profiles concurrently with the acoustic measure- 
ments. Data repeatability of ensemble-average Overall Sound Pressure Levels was 
good within a particular test day. Day to day average level variations existed. A two- 
dimensional ray tracing propagation model coupled with an empirically 
amplitude corrected predicted source noise directivity predicted the observed day to 
day average variability trends. Future research is aimed at understanding short-term, 
within a day, variability. 
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INTRODUCTION 


The Los Alamos Infrasound Program has been operating since about mid- 
1982, making routine measurements of low frequency atmospheric acoustic 
propagation. Generally, we work between 0.1 Hz to 10 Hz; however, much of 
our work is concerned with the narrower range of 0.5 to 5.0 Hz. Two 
permanent stations, St. George, UT, and Los Alamos, NM, have been 
operational since 1983, collecting data 24 hours a day. For the purposes of 
this discussion, we will concentrate on our measurements of large, high 
explosive (HE) events at ranges of 250 km to 5330 km. Because our 
equipment is well suited for mobile deployments, we can easily establish 
temporary observing sites for special events. The measurements in this 
report are from our permanent sites, as well as from various temporary sites. 
In this short report we will not give detailed data from all sites for all events; 
rather, we will present a few observations that are typical of the full data set. 

The Defense Nuclear Agency sponsors these large explosive tests as part of 
their program to study airblast effects. A wide variety of experiments are 
fielded near the explosive by numerous Department of Defense (DOD) services 
and agencies. Our measurement program is independent of this work; we use 
these tests as energetic known sources, which can be measured at large 
distances. Ammonium nitrate and fuel oil (ANFO) is the specific explosive 
used by DNA in these tests. Table I gives the test names, dates, charge 
weights, and number of our infrasonic stations operated for each test. All 
tests were fired at White Sands Missile Range, NM. 

BACKGROUND 

The basic sensor for our work is the Globe 100 microphone. A series of 
porous hoses is used to reduce the noise from the low-level local wind. 

Figure 1 shows a microphone and associated noise-reducing hoses, which can 
be thought of as a modification of the Daniel's tube used for lower frequency 
work (reference 1). During periods of quiet background, this sensor can 
easily detect signals down to a few tenths of a microbar. In our frequency 
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domain, there are natural as well as man-made sources of infrasound, some of 
which are described in Chapter 9 of reference 2. However, this background 
is not saturated with confusing signals, which simplifies the detection 
problem considerably. 

An infrasound array consists of 3 to 6 sensors placed in a regular pattern. 
We employ standard, time-delay, and sum beamforming techniques to 
process the recorded data. The present algorithm is a modified version of 
one due to Young and Hoyle (reference 3). Generally, 20-seconds of data are 
processed at a time, followed by a 50% shift and continued processing. For 
each 20-second window, the beamformer provides the correlation coefficient, 
trace velocity, and azimuth of the highest correlation signal, as well as the 
power spectrum for that interval. Longer intervals of data can be summarized 
in the manner illustrated in figure 2, where 60 minutes of data are shown. 
The presence of a signal is easily seen as the fixed azimuth line from 16:36 UT 
to 16:43 UT. 

Signal energy propagates in the atmospheric sound ducts created by the 
ambient temperature structure, or by a combination of temperature and 
wind. When propagation is in the same direction as the upper atmospheric 

winds, total refractions occur between 40 km and 60 km altitudes. The upper 

atmospheric winds are seasonal in nature, blowing to the east in winter and 

to the west in the summer (reference 4). It is important to note that at these 

altitudes the wind speed can be a significant fraction of the sound speed; 
therefore, the wind profile must be included correctly in any calculational 
work. The simplification of an effective sound speed profile is not 
appropriate for these propagation paths. 

OBSERVATIONS 


Before discussing specific time series for two events, a few general 
comments will be useful. We use the concept of an average velocity to 
broadly classify the observed propagation paths. Here, average velocity is 
just the great circle source to receiver distance divided by travel time. 

With wind propagation, the strongest signals arrive with an average velocity 
of 0.29 km/s. Ray tracing results confirm that this corresponds to total 
refractions at 40 km to 60 km altitudes. Higher average velocities indicate 
lower paths, and conversely. The stations north and northwest of White 
Sands, Los Alamos and St. George, often observe a first arrival with an 
average velocity of 0.34 km/s. This must be energy that travels at or very 
near the surface; we will for the moment refer to this as the surface wave. 
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The seasonal winds have a significant effect on observed pressures. A 
propagation path in the direction of the wind will result in a larger pressure 
than along a path directed against the wind. We apply a correction for this 
effect, which normalizes amplitudes to a zero wind condition. This procedure 
is described in reference 5. 

For the Minor Scale test we had two arrays at Barking Sands, Kauai, HI, at a 
distance of 5330 km, our most distant detection. One array operated at our 
standard frequency range, while the other operated at lower frequency, about 
0.01 Hz to 0.1 Hz. Both arrays detected the event, but the detection at the 
lower frequency was much better, as expected for such large distances. 


Figure 3 presents single channel time series from four sites for the Misty 
Picture event. Each panel is composed of 12 minutes of data in three, four 
minute windows. The time above each window is the time at the start of the 
window. Average velocities, m/s, are written below specific features for easy 
reference. For Los Alamos, note the surface wave at 342.5 m/s. For these 
time series, the data are well correlated with the source azimuth from first 
arrival to the end of the record shown. 

Figure 4 is the same as figure 2 but is for the Misers Gold event. The 
surface wave at Los Alamos is evident with a 339.2 m/s average velocity. 

Again signal energy is well-correlated with source azimuth from first arrival 
until the end of the displayed record. These two figures illustrate the 
character of the observations for these energetic events. Strong multiple 
arrivals are common, with total durations on the order of 10 minutes. The surface 
wave is common at 250 km north, and has been observed at 750 km to the 
northwest. 

In figure 5, power spectra contours are shown for the Misty Picture event 
as observed at St. George, UT. Contours of power are given as functions of 
time and frequency. Four major arrivals are seen from 18:56:30 to 19:03:00, 
with the largest powers concentrated below 1.2 Hz. Note that the major 
arrivals have frequency contributions across the whole band, from 0.2 Hz to 
3.0 Hz. 


For the purpose of examining pressure as a function of range, we have 
found it useful to place all the data on a common scale by the use of scaled 
range. In figure 6, we give the peak-to-peak amplitude of the largest signal 
(wind corrected) as a function of scaled range. The scaled range is the actual 
range divided by (2W) 1 / 2 , where W is the charge weight in tons, the factor of 
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2 is the standard factor for surface explosions, and the exponent is that 
appropriate for the cylindrical geometry of the ducted paths. The least 
squares slope is -1.4, showing only a modest increase in loss over the 
cylindrical value of -1.0. 

During the conference, a few of the participants (Drs. Raspet, Attenborough, 
and West) suggested that the surface wave was likely a creeping wave, as 
described by Pierce (reference 6). Following the discussion in reference 6, 
we have estimated the attenuation for such waves and find an attenuation 
coefficient of 0.1 km' 1 . Over the shorter path of 250 km, this gives a huge 
loss, sufficiently large, we believe, to rule out this explanation. 

We wish to acknowledge the support of the Department of Energy Office 
of Arms Control for the work supported here. 


TABLE I - EXPERIMENTS 


Event 

Date 


Sites 

Millrace 

9/16/81 

600 

1 

Pre Direct 
Course 

10/7/82 

24 

2 

Direct Course 

10/26/83 

600 

4 

Minor Scale 

6/27/85 

4800 

4 

Misty Picture 

5/14/87 

4800 

5 

Misers Gold 

6/01/89 

2400 

8 
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Figure 1. Infrasound microphone and noise reducing porous hoses. 
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Figure 2. Summary of beamformer array processing program showing 

correlation coefficient, trace velocity, and azimuth as functions of 
time. Data are for Miser's Gold event observed at Readley, Ca. 
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Figure 3. Samples of observed Misty Picture time series from four stations. 
The Bishop and Bakersfield stations were in California. 
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Figure 4. Samples of observed Miser's Gold time series from three stations. 
Crescent Valley is in Nevada, and Murphy is in Idaho. 
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Figure 5. Power spectra contours for Misty Picture from the St. George, UT 
array. 



Figure 6. Peak to peak amplitude vs scaled range for the six ME events. 
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THE CORRECTION OF INFRASOUND SIGNALS FOR UPPER ATMOSPHERIC WINDS 


J. Paul Mutschlecner and Rodney W. Whitaker 
Los Alamos National Laboratory 


INTRODUCTION 

Infrasound waves propagate in the atmosphere by a well known 
mechanism produced by refraction of the waves, return to earth, and 
reflect at the surface into the atmosphere for subsequent bounces (see e.g. 
reference 1). Figure 1 illustrates this phenomenon with results from a ray 
trace model. In this instance three rays are returned to earth from a region 
centered at about 50 kilometers in altitude and two from a region near 110 
kilometers in altitude. The control of the wave refraction is largely 
dominated by the temperature-height profile and inversions; however, a 
major influence is also produced by the atmospheric wind profile. Figure 2 
illustrates the considerable ray differences for rays moving in the wind 
direction (to the right) and in the counter direction (to the left). It obviously 
can be expected that infrasonic signal amplitudes will be greatly influenced 
by the winds in the atmosphere. The seasonal variation of the high altitude 
atmospheric winds is well documented (see e.g. reference 2). Figure 3 
illustrates this with average statistics on the observed zonal wind in the 
region of 50 ± 5 kilometers in altitude. The results are based upon a survey 
by Webb (reference 2); Webb terms this parameterization the Stratospheric 
Circulation Index (SCI). The very strong seasonal variation has the ability to 
exert a major seasonal influence on infrasonic signals. It is our purpose to 
obtain a method for the correction of this effect. 

METHODOLOGY 

There are two possible approaches to the determination of a procedure for 
the correction of infrasound signals for the effects of winds. The first of these 
is by modeling of infrasonic propagation in the presence of various wind 
profiles. We are currently taking this approach with both a ray-trace model 
and a normal mode model and hope to show results in the near future. The 
second approach is to derive a correction method empirically from a 
sufficiently large and consistent data set. This is the method which we report 
upon here. The results given here are preliminary in nature and we present 
only a simplified outline of the procedures. As indicated in Section IV, more 
extensive work in the near future will provide comprehensive results. In the 
meantime, we have been applying the results to our measurements (see 
reference 3). 
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A large data set, which appears to be appropriate for the empirical work, is 
given in reference 4. It consists of infrasonic observations by the Sandia 
National Laboratory of atmospheric nuclear tests conducted at the Nevada 
Test Site (NTS) during the period 1951 to 1962. The observations were made 
at nine stations surrounding NTS as shown in figure 4. While several of the 
stations are probably too close to the source region to be useful, at least six 
stations appear to be appropriate. A total of 80 events are presented by Reed 
and cover an explosive yield range from 1/2 to 74 kilotons (HE equivalent). 
The measurements were made by a standard set of microbarographs. 

This consistently measured and analyzed set of signals, observed at many 
times of the year, presents a unique set of data for our purposes. 

ANALYSIS AND RESULTS 

Since the atmospheric nuclear data are for a variety of yields, and also 
have effects due to some variation in height of burst, it is necessary to scale 
all data for these two factors. Reed has done this with a W^*4 scaling law (W 
= kilotons HE equivalent) and a height-of-burst functional relation. While 
some discussion of both of these scaling relations is appropriate, the 
preliminary nature of the present work leads to use of Reed’s corrections. 

Figure 5 illustrates the resulting scaled data for the St. George, Utah, 
station as a function of date. The symbols refer to various test series, 
unimportant here, and the line is an eye-fitted relation. Amplitudes are 
scaled to a 1 kiloton explosion. The very strong seasonal effect is the major 
feature of the data; the much lower amplitude during the summer period 
presumably results from the winds contrary to wave propagation to 
St. George during summer. It is unfortunate that there are no data for the 
period in mid winter, January - December. Examination of Reed's results 
shows that the seasonal variation changes markedly with station direction 
from the test site: northern and southern stations show a much weaker effect 

than does St. George. This follows from the fact that those stations are only 
slightly affected by the zonal wind and are primary affected by the much 
weaker and less variable meridional winds. This supports the hypothesis 
that we see primarily the seasonal wind effect in the data variation. 

High altitude winds are conventionally measured by rocketsonde (and now 
by satellite). The absence of rocketsonde observations for the period covered 
by the NTS data leads us to use the statistical SCI data of reference 2 as a 
first-order estimator of the atmospheric wind. The eye-fitted lines for each 
station (e.g. figure 5) were used to estimate signal amplitudes at monthly 
intervals. Figure 6 illustrates the result for St. George, where the wind 


144 


velocity component directed from NTS to St. George is plotted against the 
amplitude. Similar relations are determined for each station, using the wind 
vector toward the station based upon meridional and zonal SCI values. The 
results are reasonably consistent among the stations; however, the nearly east 
or west stations, such as St. George, are probably more dependable because of 
the strong zonal wind variation there. 

The averaged result of the analysis is, 

A = A 0 10 kv 

where A = observed amplitude (mbars) 

A 0 = corrected "zero-wind" amplitude (mb) 
k = 0.018 s/meter 

V = SCI vector from source to observer meters/s. 

The relation permits us to correct all observations for the wind to derive 
consistent "zero-wind" amplitudes. 

COMMENTARY ON RESULTS 

We have applied the method described here to a wide variety of 
observations. Where possible, we have used rocketsonde or satellite 
observations of the wind profile, deriving from these an effective SCI wind 
vector amplitude. Since the actual wind can vary widely from the statistical 
SCI, use of statistical values is less dependable. As an example of the use of 
the method, we show its application to a set of observations of signals from a 
set of earthquakes which were observed at our St. George array, with the 
exception of the earthquake at Mb = 7.8, which is from reference 5. All 
amplitudes have been scaled to a distance of 1000 kilometers by use of the 
factor (R/IOOO)!-!^ where R is the actual range in kilometers. Figure 7 
shows the amplitudes uncorrected for wind against the body magnitude, Mb, 
for each earthquake. The bars on some observations indicate the range of 
interpretation of peak amplitude. Figure 8 shows the same set of data but 
with the amplitudes corrected by our method; only statistical SCI winds were 
employed. Clearly the effect on the relationship is very large; there is also 
the impression that the connection between amplitude and Mb may be 
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clearer. However, a mechanism is not now known, so this evidence is 
circumstantial. The smallest earthquake, at Mb = 4.6, is an aftershock and 
may not be in the same class as the others. 

We have not completed our analysis of the nature of the wind effect on 
amplitudes and therefore can only speculate on the correction method and its 
form. It is clear that two mechanisms can change the amplitude with change 
in the wind profile: (1) the number of rays captured at a return layer will 
vary depending upon the wind, (2) the distribution of rays on the surface will 
depend on the wind profile (e.g., see figure 2). Modeling should help in the 
understanding of the correction form. 

It should be noted that all of the work reported here applies to those 
signals resulting from returns from a 50 kilometer high level. We believe 
that signals resulting from the 110 kilometer level require an additional 
correction and have formulated a tentative correction formulation. Such 
signals, in general, will occur only with near zero wind, or counterwind, 
conditions. 


FUTURE EFFORTS 

We have reported here on our preliminary results. We are now working 
on an improved analysis. This will include the following: 

1 . A comprehensive statistical analysis of the data set resulting in an 
improved formulation; 

2. Modeling of the effect to better understand its physical basis; 

3. A detailed investigation of the "counterwind" signal circumstances, 
using appropriate data and modeling. 


We wish to acknowledge the support of the Department of Energy Office of 
Arms Control for the work supported here. 
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Figure 1. Results of a ray trace model showing five sample rays. Rays 1, 2, 3 
are refracted and returned from a region near 50 kilometers and 
rays 4 and 5 from near 110 kilometers. 
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Figure 2. Ray tracing results for propagation with the wind (to the right) and 
against the wind (to the left). 
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Figure 3. This illustrates the seasonal variation of the average SCI wind at 
White Sands. 
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Figure 4. 
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Figure 5. Nuclear atmospheric burst scaled signal amplitudes recorded at the 
St. George, Utah, microbarograph station against date. 



Figure 6. Monthly average values of nuclear atmospheric burst scaled signal 
amplitude versus statistical SCI velocity for St. George. 
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I. INTRODUCTION 

Measurements of sound propagation over distances up to 1000 m were carried out with an impulse 
sound source offering reproducible , short time signals. Temperature and wind speed at several heights 
were monitored simultaneously; the meteorological data are used to determine the sound speed gradients 
according to the Monin-Obukhov similarity theory. The sound speed profile is compared to a 
corresponding prediction, gained through the measured travel time difference between direct and 
ground reflected pulse (which depends on the sound speed gradient). Positive sound speed gradients 
cause bending of the sound rays towards the ground yielding enhanced sound pressure levels. The 
measured meteorological effects on sound propagation are discussed and illustrated by ray tracing 
methods. 


II. SOUND SPEED PROFILES 

Wind speed and temperature are functions of elevation above ground. They are interrelated and can be 
described by the Monin-Obukhov similarity theory /l/ using two parameters, the friction velocity u* 
and the scaling temperature T*. The Monin-Obukhov Length L is a stability parameter for the turbulent 
atmospheric surface layer 

L = (T m /gkXu* 2 /T*) 

T m : representative temperature, g: acceleration due to gravity, k: von Karmans constant (0.41). 

The sound speed profile can be described by 

j 

c(z) = c(z 0 ) + a’(/n(z/z 0 ) + V<z/L)) (1) 

| where a’= u*/k + 0.6T*/k and z Q is the roughness length. 

1 For the stable case (positive temperature gradient, strongest bending of sound rays towards the ground), 
V<z/L) = 5z/L, and the sound speed gradient is: 

dc/dz = (a’/z)(l+5z/L) (2) 

Examples of measured temperature and wind speed profiles are shown in Fig. 1. The parameters u* and 
T* are calculated from those measured data by least square methods /2/. In the unstable case wind 
speed and temperature almost remain constant in larger elevations (no sound speed gradient), while in 
the stable case there is still an increase in wind speed and temperature. Close to the ground the profiles 
are ’logarithmic’ in both cases. 
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III. ACOUSTICALLY MEASURED SOUND SPEED GRADIENTS 

Measured travel time differences between direct and ground reflected sound can be used to estimate the 
sound speed profile during the measurement. The geometrical time difference At g is: 

At g = AD/c = [/(h g + h r ) 2 + D 2 ’ - v/(h g - h r ) 2 + D 2 "]/c (3) 

c: sound speed ; h g : source height; h : receiver height. 

The sound speed profile c(z) causes an additional time difference At’, the total time difference At 
between direct and ground reflected sound is: 


At = At’ + At (4) 

S 

For small deviations of the actually curved ray path from the geometrical path the travel time of the 
ground reflected pulse can be estimated to be for the stable case /3/: 

T refi = J c ( z )’ lds = Dc(z o r 2 [c(z 0 ) -a7«(z’ /z 0 ) + a’ - 2.5a’z’/L] + At g (5) 

where z’ = h g = h r . The travel time of the direct pulse at height z’ is: 

T d . r = D/ c(z’) = D/[c(z 0 ) + a7«(z’/z 0 ) + 5a’z’/L] 

= Dc(z q )' 2 [c(z 0 ) - a/n(z’/z 0 ) - 5a’z’/L] (6) 

and the travel time difference is: 

At = T refl - T d . r = [a’(l + 2.5z’/L)]D/c(z 0 ) 2 + At g . (7) 

Two measurements of the travel time difference at different heights z’ are necessary to determine the 
parameters a’ and L. 

Close to the ground or for the nearly neutral case (z/L becomes 0) the sound speed profile becomes 
’logarithmic’, 

c(z) = c(z 0 ) + a/n(z/z 0 ) and dc/dz = a/z, (8) 

the parameter a can be calculated from the measured time difference for only one source-receiver 
height: 

a = At’c(z Q ) 2 /D (9) 
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IV. SOUND PROPAGATION MEASUREMENTS 

The impulse sound source /4/ used for the acoustical measurements consists of a capacitor (100 /xF, 
3.5 kV), which is discharged over a spark gap. A reproducible short pulse (less than 1 ms) is radiated 
spherically (sound pressure level 150 dB at 1 m distance). As an example Fig. 2 shows a measured time 
signal close to the sound source. The time delay between direct and ground reflected pulse (here 8 ms) 
is mainly caused by geometry. A reference signal monitored in an unechoic room at 6.25 m distance to 
the sound source is used to calculate the SPL re free field for each frequency. 

Two examples of downwind sound propagation measurements (8 August 1988 and 3 November 1988) are 
shown in the following. Measurements for several geometries done in the afternoon (unstable conditions) 
were repeated a few hours later under stable conditions (positive temperature gradient). The resulting 
meteorological effects are discussed. 

August-measurement 

For a distance of 250 m the measured time signal is shown in Fig. 3 a (temperature and wind speed 
profile see Fig. 1 c). The travel time difference At resulting from geometry is 3.7 ms. The measured 
time difference At (determined from the magnitude of the autocorrelation function of the measured 
sound pressure signal) is 5.3 ms (see Fig. 3b), the additional time delay At’ = 1.6 ms being due to 
meteorological effects (wind speed and temperature). It can be used (eq. 8,9) to calculate the sound 
speed gradient: 

a = At’c(z 0 ) 2 /D= 0.7 m/s; dc/dz = a/z = 0.7 1/s 

A few hours earlier (meteorological conditions see Fig. 1 b) a time difference At = 4.6 ms was measured 
for the same geometry, yielding a parameter a = 0.4 m/s. 

For a source and receiver height of 12.5 m there is only a small change in the measured SPL due to 
meteorological effects (Fig. 4 a). The interference pattern is shifted a little towards lower frequencies 
for the evening measurement, the sound pressure level increasing about 1 dB. For source and receiver 
situated closer to the ground (2 m. Fig. 4 b) the different meteorological conditions yield an evident 
shift of the ’ground dip’ to lower frequencies. For a distance of 1000 m (h : 12.5 m and 5 m; h : 5 m 

. 8 ’ r ’ 

F i g. 5) no ’ground dip’ occurs in the measured frequency range. 

N ovember-measurement 

Fig. 6 a-d show SPL’s for different geometries and two meteorological conditions. For h g = h r = 5 m 
and 100 m distance (SPL see fig. 6 a) a time difference At = 1.7 ms was measured at 13.40, increasing to 
At = 2 ms two hours later (wind speed about 2 m/s in both cases, but negative temperature gradient in 
the first and positive temperature gradient in the second case). The time difference At g due to geometry 
is 1.5 ms. At’ increases from 0.2 ms to 0.5 ms, the sound speed gradient in the late afternoon is more 
than twice as large as in the early afternoon. 

The meteorological effects are strongest for low source heights and for large distances. The ’ground dip’ 
is shifted towards lower frequencies with increasing sound speed gradient. If source and receiver are 
closer to the ground than in Fig. 6 a (h g = h. = 1.5 m, D = 100 m. Fig 6 b) an evident ’ground dip’ 
around 400 Hz occurs which is reduced for the stable measurement. Fig. 6 c shows a measurement at 
825 m distance (h g = 1.5 m, h = 5 m). The broken line is the calculated SPL using a single- parameter 
impedance model (see /5/) and a flow resistivity typical for grass covered ground (no sound speed 
gradient assumed). The positive sound speed gradient reduces the ground dip. For the stable case 


157 



(strongest sound speed gradient) the SPL is increased some 20 dB around 500 Hz. Fig. 6 d shows the 
SPL for a source height of 5 m. The measured SPL is larger than that for a source height of 1.5 m. 

A ray tracing simulation for the 825 m measurement under stable conditions (Fig. 6 d) is shown in 
Fig. 7. Multiple reflections at the ground and focussing effects (see /6/) should result in larger sound 
pressure levels. Indeed, the sound pressure level is increased about 10 dB for the stronger sound speed 
gradient. The measured time signals corresponding to the spectra in Fig. 6 d also show the expected 
differences (see Fig. 8). While for the ’unstable’ measurement direct and ground reflected sound arrive 
almost at the same time, a lot of ground reflected pulses with enhanced pressure levels are observed for 
the ’stable’ measurement yielding a signal of more than 10 ms length. Calculations with ray tracing 
methods /7/ for the stable meteorological condition predict a travel time difference between the direct 
ray (which arrives first) and the latest multiply reflected ray of 13 ms, in good agreement with the 
measurement. 

In Fig. 9 the sound speed gradient is plotted as a function of elevation above ground. Curve (a) 
represents the gradient calculated from the acoustical measurement, curve (b) is a best fit to the 
measured wind speed and temperature values. Good agreement is achieved close to the ground where the 
sound speed profile has a ’logarithmic’ shape. With increasing height the curves differ due to the 
influence of the stability (z/L is not small compared to 1 in eq. (2)). 


V. CONCLUDING REMARKS 

Sound speed gradients determined from acoustical measurements represent integrated values during the 
actual travel time and along the actual sound path. Meteorological data in order to determine sound 
speed gradients, on the other hand, are measured locally and require several sensors (wind speed and 
temperature) of high accuracy close to the propagation area. 

Downwind sound propagation (comparable wind speeds) is extremely sensitive to the stability of the 
atmospheric surface layer. Positive temperature gradients (stable conditions) yield a positive sound speed 
gradient even at large elevations, where negative temperature gradients (unstable conditions) yield a 
negligible sound speed gradient. Close to the ground the sound speed profile has a ’logarithmic’ shape 
and the sound speed gradient can be described by one parameter a, which can be calculated from the 
measured travel time difference between direct and ground reflected sound. 

The bending of sound rays towards the ground is strongest under stable conditions. The ’ground dip’ is 
diminished and shifted to lower frequencies yielding negligible excess attenuation in the frequency 
range relevant for noise propagation. Focussing effects and multiple reflections lead to enhanced sound 
pressure levels. 
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Li Measured wind speed and temperature profiles in the atmospheric surface layer (+). The solid 
line fits the measured values using the Monin-Obukhov similarity functions. 
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Fig. 3 a: Measured time signal 
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Fig. 3 b: Autocorrelation function 
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Fig. 9: Sound speed gradient [ 1/s] 
as a function of elevation above ground 

a) measured acoustically (Fig. 3) 

b ) determined by meteorological measurements 

(Fig. 1 c) 
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Abstract 

The effect of refraction due to wind and temperature gradients on energy received 
from low flying aircraft is examined. A series of helicopter and jet flyby's were recorded 
with a microphone array on two separate days, each with distinctly different meteorological 
conditions. Energy in the 100-200 Hertz band is shown as a function of aircraft range 
from the array, and compared with the output of the Fast Field Program. 


I. Introduction 

This paper examines the effect of wind and temperature gradients on energy received 
at a microphone array from a series of aircraft flyby's. Of interest is the energy contained 
between 100 and 200 Hertz, the frequency band used in our acoustic detection and tracking 
algorithms. 

One aspect of this work is to estimate our ability to detect and track low flying 
aircraft, or conversely, to assess the vulnerability of aircraft to acoustic detection and 
tracking. Propagation characteristics, which are largely influenced by wind and 
temperature gradients, must be taken into account if we are to make accurate predictions. 

To illustrate the impact that wind and temperature gradients can have, received energy 
as a function of aircraft range has been calculated from aircraft flyby's on two separate 
days, each with distincdy different meteorological conditions. Sound speed profiles, 
derived from wind and temperature data collected during the experiments, are used to 
generate ray plots. Visualization of the ray paths helps to explain features seen in the 
experimental data. 

To predict detection range or tracking ability for a given set of meteorological 
parameters, we must estimate acoustic energy as a function of distance from the source. To 
this end, the output of a propagation model, the Fast Field Program, is compared to the 
experimental results. 


* This work was sponsored by the Department of the Air Force. 
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II. Experiment 


Aircraft flyby's, depicted in Figure 1, were recorded on two different days 
(designated as Day 1 and Day 2). Results presented here are from a helicopter on Day 1, 
and a jet aircraft on Day 2. Both aircraft flew in a straight line at a constant altitude past a 
nine element microphone aiTay. Ground truth TSPI (Time SPace Information) of the 
aircraft's position and velocity, corrected for acoustic propagation time, was also recorded 
during each flyby. Details are given in Table 1. 

Array data were sampled at 2048 samples/second during the experiment and recorded 
directly to magnetic tape. The array consisted of nine GenRad 1962-P42 microphones with 
standard Sennheiser windscreens. Microphones were placed in notched wooden blocks on 
the ground in a tri-delta configuration (reference 1). The array was used with a wideband 
direction finding program (reference 2) to aid in determining whether received energy was 
signal from the aircraft, or noise. This is discussed further in Section IV. 

Meteorological data (temperature, wind speed and direction, and relative humidity) 
were recorded to a height of 300 meters before and after the experiment using a tethered 
balloon. These parameters were also recorded on the ground throughout both 
experiments. Meteorological data were stored every 10 seconds during the experiment. 

The wind was from the South (190 degrees) on Day 1, and from the North (15 degrees) on 
Day 2. Headings of 345 and 165 degrees put the aircraft approximately into the wind, or 
with the wind. 

The helicopter was louder when it was inbound to the array, so only incoming 
portions of the helicopter data are analyzed. There were two runs incoming from the North 
(345 degrees), and two runs incoming from the South (165 degrees). The closest point of 
approach (CPA) from each direction was 90 and 230 meters. 

The jet was louder outbound from the array, so only outgoing portions of those runs 
are used. Three runs outgoing to the North (345 degrees), and three runs outgoing to the 
South (165 degrees) are analyzed. The CPA for these runs varied from 140 meters to 716 
meters. 


III. Data 

Array data 

The array time series for one of the helicopter runs at its CPA is shown in Figure 2a. 
This same time series is displayed in Figure 2b after bandpass filtering between 100 and 
200 Hertz. The spectra from two of the channels are shown in Figure 3. These spectra 
show the strong harmonic structure that is typical for helicopters. 
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Array time series for the jet are shown in Figures 4a and 4b. The jet spectrum from 
two of the channels are shown in Figure 5. These figures illustrate the broadband spectra 
that is typical of jets. 

The drop in power level in both Figures 3 and 5 at about 750 Hertz is due to the 
antialiasing filter. A rise in energy below 50 Hertz in the spectra of Figure 5 is due 
to wind noise. 

Environmental data 

Meteorological data collected from a tethersonde was used to calculate sound speed as 
a function of height Data taken during one of the balloon raisings on Day 1 is shown in 
Figure 6. There was a normal temperature lapse above 50 meters, with the wind out of the 
South. Sound speed profiles for Day 1 at 345 degrees (looking North of the array) and 
165 degrees (looking South of the array) are shown in Figure 7. 

On Day 2, the wind was from the North (Figure 8). The wind speed initially 
increased up to 70 meters, then decreased with height, up to 300 meters. This unusual 
wind profile, along with a temperature inversion, led to the sound speed profile in figure 9. 


IV. Analysis 

Energy as a function of range 

Received energy is calculated for each one-second segment (2048 points) of the array 
time series. This corresponds to a spatial average of about 30 meters for the helicopter, and 
240 meters for the jet The power spectrum for each channel is first calculated using a 
Hanning window and 2048 point fft's. After integrating the power spectra between 100 
and 200 Hertz, the values for all channels are averaged. The level calculated for that one- 
second segment is then matched to the corresponding TSPI range, yielding energy received 
at the array when the aircraft was at that particular range. 

Separating signal from noise 

It is not always clear if acoustic energy received at a microphone is signal from an 
aircraft, or wind noise. Whether it is signal or noise will depend upon the propagation 
conditions (for example, the presence of a shadow zone), the level of wind noise, and the 
distance from the aircraft to the microphone. Discriminating between signal and noise is 
important when comparing the output of a propagation model to experimental data; we do 
not want to ascribe propagation effects to our experimental data when no signal is there to 
model. To ensure that we were only looking at signal from the aircraft, the array time 
series was used with a wideband direction finding algorithm (reference 2) to classify the 
received energy as signal or noise. 

The direction finding algorithm outputs the energy arriving along a specified number 
of directions. The direction from which the maximum energy arrives is the detected azimuth 
of the source. Energy and azimuth pairs for other directions are output in order of 
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decreasing received energy. For energy to be selected as signal from the aircraft, we 
require the detected azimuth to be close to the azimuth reported by the TSPI (ground truth) 
data. In addition, we require that energy coming from the direction of the detection be 
larger than energy coming from other directions, otherwise we are probably measuring 
ambient noise. All energy versus range data reported in the next section have been 
screened using the above criteria. 

Received energy data 

To help in understanding features in the received acoustic energy data, raytraces were 
calculated (reference 3) using the sound speed profiles from Figures 7 and 9, and are 
shown along with the energy versus range plots. The ray plot for the case when the 
helicopter was incoming from the North on Day 1 (calculated from the sound speed profile 
in Figure 7a) is shown in Figure 10a. 

If the aircraft is considered to be at zero range and an altitude of 40 meters on the ray 
plot, then the number of rays intersecting the ground at any range gives an indication of the 
acoustic energy that would be received at that distance from the aircraft. Since the sound 
speed decreases with height (Figure 7a), rays leaving the aircraft bend upward, and a 
shadow zone is formed at about one kilometer from the source. 

The received acoustic energy as a function of range for runs in which the helicopter 
was incoming from the North is given in Figure 10b. Each data point represents the energy 
averaged over one second in the 100 - 200 Hz. band. To provide a reference, a solid curve 
representing spherical spreading is shown along with the experimental data. As suggested 
by the raytrace, there is a larger decrease in received energy than predicted by spherical 
spreading past one kilometer, where the shadow zone begins. Note that the energy level 
drops significantly in the shadow zone, but does not go to zero, as ray theory predicts. 

The raytrace and energy plot for runs in which the helicopter was incoming from the 
South are shown in Figure 1 1. In this case, the sound speed increased with height (Figure 
7b), causing the rays to be bent downward. Energy received past about one kilometer is 
less than that predicted by spherical spreading since much of the energy is refracted 
downward at short ranges; rays are more spread out at longer ranges than would be the 
case for spherical spreading. Other factors, such as directivity of the source, and the 
ground effect, are likely to be present as well. 

The raytrace (calculated from the velocity profile in Figure 9a) and energy plot for 
outgoing runs to the North on Day 2 (jet) are given in Figure 12. The raytrace suggests a 
received energy somewhat higher than indicated by spherical spreading at short ranges 
where the rays are refracted downward, and less received energy at longer ranges where 
the rays are refracted upward. Comparison of the experimental data and the spherical 
spreading curve shows this to be the case. 

When the aircraft was South of the array, an initial decrease in sound speed up to 80 
meters in height (Figure 9b) caused shallow angle downgoing rays to be bent upward, 
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creating a small shadow zone. Past 80 meters, there was a general increase in sound speed 
with increasing height, which caused the rays to be bent downward. The steep drop in 
received energy (Figure 13b) between one and three kilometers corresponds to the shadow 
zone seen in the raytrace. There is an increase in energy between four and six kilometers as 
rays leaving the source with an upward angle were refracted back downward. 


V. Comparison with FFP 

Sound speed profiles in Figures 7 and 9 were used as input to the Fast Field Program 
(references 4-6). As seen in Figures 14 and 15, agreement between the model output and 
general features in the experimental data is quite good. In particular, note that the FFP 
output closely models the experimental data in the shadow zones seen beyond one 
kilometer in Figure 14a, and between one and three kilometers in Figure 14b. 


VI. Summary 

Measurements of acoustic energy from a series of aircraft flyby's were presented. 
Features in the experimental data were explained in terms of the propagation characteristics 
present at the time. Sound speed profiles, from meteorological data taken during the 
experiment, were used as input to the Fast Field Program. The FFP was seen to provide 
an excellent prediction of the general features found in the experimental data. The large 
difference between the experimental results and simple spherical spreading emphasizes the 
need for accurate and detailed meteorological data. 
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Table 1. 


Day 

Aircraft 

Heading 

(degrees) 

Velocity 

(m/sec.) 

Wind speed 
(m/sec.) 

Wind direction 
(degrees) 

1 

helicopter 

345/165 

30 

1.3 

190 

2 

jet 

345/165 

240 

2.9 

15 



Figure 1. Aircraft flyby. 
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Figure 2. Helicopter time series at closest point of approach (90 meters). 


172 


SPL (dB re 20 uPa) SPL (dB re 20 jiPa) 


Channel 2 



Channel 8 



Figure 3. Helicopter spectra at closest point of approach (90 meters). 
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Figure 4. Jet time series at closest point of approach (716 meters). 
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Figure 5. Jet spectra at closest point of approach (716 meters). 
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a) North of the array (345 degrees). 
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b) South of the array (1 65 degrees). 


Figure 7. Sound speed profiles on Day 1 . 
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Figure 9. Sound speed profiles on Day 2. 
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a) raytrace calculated from the velocity profile in Figure 7a. 
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Figure 10. Helicopter incoming from the North (Day 1). 
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a) raytrace calculated from the velocity profile in Figure 7b. 



b) experimental data. 


Figure 1 1 . Helicopter incoming from the South (Day 1). 








Energy (dB re 20 (iPa) Height (meters) 


Range (km) 

01 23456789 

400i * * . * — — , 



a) raytrace calculated from the velocity profile in Figure 9b. 



b) experimental data. 


Figure 13. Jet outgoing to the South (Day 2). 
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a) helicopter incoming from the North. 



b) helicopter incoming from the South. 


Figure 14. Comparison of Day 1 experimental data with the FFP. 
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b) jet outgoing to the South. 


Figure 15. Comparison of Day 2 experimental data with the FFP. 
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COMPARISON OF FFP PREDICTIONS WITH MEASUREMENTS OF A 
LOW-FREQUENCY SIGNAL PROPAGATED IN THE ATMOSPHERE 

D. Keith Wilson and Dennis W. Thomson 
Department of Meteorology, 

The Pennsylvania State University 

SUMMARY 

An experimental study of low-frequency propagation over a distance of 770 m was 
previously reported [J. Acoust. Soc. Am. Suppl. 1 86, S120 (1989)]. For that study, 
sound speed profiles were reconstructed entirely from surface-layer micrometeorological 
data. When the acoustic data were compared with theoretical predictions from a fast field 
program (FFP), it was found that the FFP underpredicted sound levels measured in a 
shadow zone. In this paper, the effect on the predictions of including meteorological data 
for heights greater than the surface layer, i.e., wind profiles measured by a Doppler sodar, 
is discussed. Vertical structure of turbulence is simulated by stochastically perturbing 
the mean profiles, and the agreement between the acoustic data and FFP predictions is 
improved. 


INTRODUCTION 

Previous studies of fluctuations in acoustic signals propagated in the atmosphere have 
typically been concerned with time scales of a few seconds or less. The purpose of such 
experiments was primarily to study scattering by turbulence with sizes on the order of the 
acoustic wavelength. 

The experiment described in this paper was designed to study temporal variability on 
a much longer time scale. The level of a low-frequency signal was monitored for several 
periods lasting between two and six days, with the sound level being recorded at one 
minute intervals. The atmospheric phenomena affecting acoustic signals on these time 
scales are large-scale turbulence (e.g., thermals), diurnal evolution of the atmospheric 
boundary layer, and synoptic-scale weather systems. 

Along with monitoring of the acoustic signal, a wide variety of micrometeorological 
data were logged. One use of these micrometeorological data was the reconstruction of 
half-hour mean sound speed profiles. The sound speed profiles were used in a propagation 
model, the fast field program (FFP). 

In an earlier paper presented at the fall 1989 meeting of the Acoustical Society of 
America [1], comparisons of acoustic data with predictions from the FFP were presented. 
The agreement between the data and predictions was found to be reasonably good, so long 
as the receiver was not in a shadow zone. For a receiver in a shadow, the disagreement 
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was up to 20 dB. 

In the first section of this paper, the earlier paper will be summarized. In particular, the 
experimental procedures and the method originally used in the profile reconstructions will 
be discussed. In the second section, the profile reconstruction is extended to include data 
recorded by a Doppler sodar. The new method uses the generalized inverse to construct a 
least squares fit to the meteorological data. In the third section, stochastic perturbations 
are added to the profiles, in order to model the vertical structure of turbulence. 

SUMMARY OF PREVIOUS RESULTS 
Experimental Procedures 

All of the experimental data were collected at the Pennsylvania State University’s Rock 
Springs Agronomy Research Center. The propagation path was over crop area which had 
been planted with corn and soybeans, although the crops had been harvested prior to the 
experimental runs. 

The basic experimental plan was quite simple: a 27.7 Hz source was continuously 
monitored at a distance of 770 m by one or two microphones, which were connected to 
a computerized data logging system. The source consists of four identical boxes, each 
approximately 1 cubic meter in size and having two fifteen-inch diameter moving-coil 
loudspeakers. The boxes were made to resonate at low frequency by drilling suitably-sized 
ports. The final operating frequency of 27.7 Hz was chosen because it was approximately 
the mean of the resonance frequencies of the individual boxes. The loudspeakers were 
driven by two Altec 400-watt stereo amplifiers. The sound pressure level at a reference 
microphone (General Radio 1560), located 8 meters from the source, was measured using 
a Hewlett-Packard 3561 single-channel spectrum analyzer, and logged to floppy disk via 
an HP 9186 microcomputer. 

The set-up at the remote receiving station was quite similar. The sound pressure 
level at the two remote microphones (General Radio 1560) was monitored with a Hewlett- 
Packard 3562A dual-channel spectrum analyzer. The microphones were 0. 6-2.0 m from 
the ground. Fifteen FFT’s were performed on the microphone signals each minute and 
averaged by the spectrum analyzer. The total power in the band 27.7 ± 1.0 Hz was 
calculated with an HP 9186 microcomputer and logged to floppy disk. 

Preliminary investigations were performed in September 1988. The first extended 
experimental run took place during 13-16 October 1988. At this time the ground had not 
yet frozen, and was dry at the surface. The experiment then was repeated three times in 
February 1989, and once in March 1989. 

To facilitate comparison of the various experimental runs, the sound level at the 8 m 
reference microphone was used to normalize the data. The first step in the normalization 
process was to compute an estimated sound pressure level at a distance of 1 m from the 
middle of the source. By assuming that the spreading from the source to the reference 
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microphone was approximately spherical, the 1 m reference SPL could be found by adding 
18 dB to the 8 m measurement. All of the acoustical data at the 770 m microphone were 
then normalized by subtracting the 1 m reference level. No additional compensation was 
made for cylindrical or spherical spreading losses between the 1 m reference and the remote 
microphones. 


Meteorological Measurements 

The surface layer measurements regularly logged at Rock Springs are extensive. Only 
those measurements which are most useful for interpretation of the acoustical data are 
discussed in this section. 

One of the most sensitive and versatile instruments at Rock Springs is the Kaijo Denki 
DAT-300 ultrasonic anemometer-thermometer. This device, positioned ten meters above 
the ground, samples the wind velocity and temperature at a frequency of 20 Hz. Data from 
the ultrasonic anemometer-thermometer are also used to compute a number of turbulence 
statistics; among these are the covariance of the vertical wind component w with the 
horizontal wind component u, and the covariance of the vertical wind component w with 
the temperature T. These covariances are basic to the study of momentum and heat 
transfer processes in the surface layer. 

In addition to the ultrasonic anemometer-thermometer, many other anemometers and 
thermometers with slower sampling rates are maintained. Cup anemometers and vanes are 
positioned at 2 m and 6.4 m. A device that is particularly useful in the reconstruction 
of temperature profiles is the “temperature difference probe,” which continually senses 
the temperature difference between thermistors placed at 1.9 and 8.9 m. The thermistors 
are coupled in a bridge which ensures accurate evaluation of the temperature difference. 
This device is preferrable to the use of separate thermometers placed at different heights, 
because the latter method is sensitive to errors in the absolute calibration of the separate 
thermometers . 

The data from the anemometers and thermometers are averaged for a half hour before 
being logged. Thus, with the current procedures, a half hour is the minimum interval for 
reconstruction of sound speed profiles. 

A Doppler acoustic sounder (sodar) is also available for mapping temperature structure 
and wind profiles at heights greater than the instrumentated towers. The temperature 
structure information can be used to monitor the inversion height, to observe the presence 
of stable layers in the boundary layer, and to observe the passage of thermal plumes. The 
sodax was programmed to have a height resolution of fifty meters, with wind data being 
recorded every ten minutes. 

Reconstruction of Sound Speed Profiles 

Height-dependent sound speed profiles axe required as input to refraction models such 
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as the FFP. It is highly desirable to be able to reconstruct useful profiles from only surface- 
layer data and remote measurements. Balloon launches are comparatively expensive. 

In the literature on atmospheric acoustics, sound speed profiles have typically been 
determined in two ways. The first method, which might be called the “direct” method, 
consists of measuring the wind velocity and temperature at a large number of heights. In 
practice the direct method is, of course, costly if separate sensors are used at each height. 
It is also extremely sensitive to calibration of the individual sensors. A related option 
would be to use one set of moving sensors, although this procedure is complicated by the 
presence of turbulent fluctuations in the fields [2], 

A second and more commonly used method, which might be called the “logarithmic” 
method, consists of measuring wind velocity and temperature at two heights. The mea- 
surements are then fit with a logarithmic sound speed profile. The problem with the 
logarithmic method is its accuracy: the sound speed profile is only approximately loga- 
rithmic, unless conditions are near neutral. 

The method described here to reconstruct sound speed profiles is based on surface 
layer similarity scaling theory. It is similar to the logarithmic method, in that sensors are 
required at only two heights. With surface-layer similarity scaling, however, the analysis 
does not need to be limited to neutral conditions. One of the best summaries of surface- 
layer scaling theories is Stull [3], Chapter 9. Much of the following material is presented 
in more detail by Stull. 

The type of scaling used in this paper is due to Monin and Obukhov. The meteoro- 
logical profiles are written as functions of z/L , where z is the height from the surface and 
L is called the Monin-Obukhov length, which can be written 


L = 


g*T m ' 


( 1 ) 


In the above, g is gravitational acceleration, k = 0.4 is the von Karman constant, T 0 is 
the mean surface temperature in Kelvin, 

u* = \J —w’u’ (2) 


is the friction velocity, and 


T„ = 


w’V 


u t 


( 3 ) 


is the surface-layer temperature scale. The covariances are evaluated in the surface layer. 

Note that if w'T' is positive, which is the case for statically unstable conditions, then 
T» < 0 and L < 0. In fact, for unstable conditions, — 0.5L is approximately the height 
at which buoyant production of turbulence begins to dominate over mechanical (shear) 
production (Stull, p. 182). The limit z/L —+ 0 represents neutral conditions. When L > 0, 
conditions are stable. In this case z/L indicates the extent to which mechanical turbulence 
is suppressed by the static stability in the mean temperature profile. 
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The rates of change of mean wind and temperature are written 

dujz) u* ( z \ 

dz kz M vL/ ’ 


( 4 ) 


(5) 


where the scaling functions <f> are determined through a combination of theory and data 
regression. 

When Eqs. 4 and 5 are integrated, the results can be written: 


uW “ 7 [ ln (i) - * M (I)] ' < 6 > 

T(z) = T(zt) - lAz - *.) + £ [in (£) - <h, (|)] , (7) 

where z a is the aerodynamic roughness length , and z t is the thermal roughness length. For 
unstable conditions, L < 0, 

* u (z) = (z) /°- 74 = ( 4 - 13 z) ' ■ (8) 

When these functions are integrated, it can be shown that 


(z) = * H (f) /0 ' 74 = i ln 


(1 + ) 

3 




(9) 


For stable conditions (L > 0), the (^-functions are much simpler: 


*"(z) = 1 + 47 (z)' *"(z) =a 74 + 4 ' 7 (i)' ( 10 > 


Integration yields: 

*" (z) = ' / ’ H (z) = ~ 47 (z) ' (”> 

When a plant canopy or other roughness elements are present, the effective ground 
plane should be shifted upward by an amount d , called the displacement height. This is 
an extrapolated height at which the wind speed is approximately zero. It is possible to 
estimate z a and d from the height of the roughness elements (plant canopy, buildings, etc.). 
According to Panofsky and Dutton [4], 



x canopy height, 
x canopy height. 


( 12 ) 

(13) 
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There are now three unknowns in the wind and temperature profile equations: it*, T,, 
and the surface temperature, T(z t ). This means that three independent measurements are 
required. In the original reconstruction procedure, these consisted of the temperature at 
two heights, say z\ and z 2 , and the wind velocity at one height, z 3 . Solving Eq. 6 for it, 
and evaluating at z 3 , we have 


Kujza) 

ln[(z 3 - d)/z 0 ] - - d)/L] ’ 


(14) 


Solving 7 for T», evaluating at z 2 and zi, and subtracting the z 2 equation from the z\ 
equation, z t is eliminated and we have 


k[AT + Jd(z2 ~ -gi)] 

ln[(z 2 - d)/(zi - d)] - %I>h[(z 2 - d)/L\ + V>//[(x i - d)/L]' 


(15) 


Since L is a function of it* and T,, these equations actually contain it* and T* on both sides. 
However, the ip functions are small compared to the logarithmic term, and the equations 
may be solved by first neglecting the ip functions, and then iterating until values for it, 
and T, have been converged upon. The solution is normally well behaved so long as there 
is good mixing in the surface layer. When local values of L axe less than about 5, this 
method does not converge to a solution. 

Once the wind and temperature profiles have been determined, an effective sound 
speed can be computed by adding the component of the wind velocity in the direction of 
propagation to the actual sound speed. An example of the sound speed profile evolution 
for 8 March 1989 is shown as Fig. 1. For this figure, L was assigned a value of 5 if there 
was no convergence to a solution. This was necessary during most of the nighttime hours. 


Comparison with Acoustic Data 

When the data shown in Fig. 1 were used in a fast field program (written by one of 
the authors, D. K. Wilson), the predictions shown in Fig. 2 result. For this figure, the 
acoustic data are plotted as half-hour means, so that the averaging periods of the acoustic 
data and meteorological measurements coincide. The sound speed profile was partitioned 
at 0.8, 1.4, 3, 6, 12 and 24 m. The ground was modelled as a rigid porous medium, with 
static flow resistivity of 200 000 mks rayls/m, tortuosity 2.5, and porosity 0.3. 

Notice that, prior to 0800 and after 2000, the FFP predictions agree fairly well with 
the experimental data. During the daytime hours, however, the agreement is very poor. 
Examination of the profile reconstructions shows that the predictions are in good agree- 
ment so long as the sound speed increases with height, i.e., the receiver is in a surface 
sound channel. When the receiver is in an acoustic shadow zone, the FFP predicts a much 
greater propagation loss than actually occurs. 
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INCLUSION OF SODAR WIND PROFILES 


The original method of reconstructing the profiles, discussed in the previous section, 
used only surface-layer data. This raises the question of whether incorporation of data for 
heights above the surface layer could improve the predictions. 

Recall that a Doppler sodar was in operation during the experiment. The sodar moni- 
tors wind profiles at heights above 50 m. In this section, a new method for reconstructing 
the profiles, which incorporates the sodar data, is discussed. Unfortunately, the recon- 
struction procedure is much more difficult when the sodar data are included. The main 
reason is that the problem is now over-determined: there are more experimental data than 
parameters in the model. 

Let us formulate the problem as follows. We arrange the meteorological data (surface- 
layer winds and temperatures, and Doppler sodar profiles), as a column vector d. The 
model parameters, tz„, T* and the surface temperature, comprise the column vector m. 
Retaining just the first term in the Taylor series, the forward problem can be written 

d' « Gm', (16) 

where the primes indicate the fluctuation about the actual value, and 

r - 

,J dm / 

The derivates can be evaluated, for example, by numerically differentiating Eqs. 6 and 7. 

What we desire, however, is a solution to the inverse problem; that is, we want the 
model parameters in terms of the data. Some operator (G) -1 must be constructed, so 
that 

m' » (G) _1 d'. (18) 

Construction of such operators belongs to the field of inverse theory. An excellent intro- 
duction can be found in Chapter 12 of Aki and Richards [5]. The problem is not trivial, 
because G is not normally an invertable matrix. 

For the profile reconstructions, a weighted generalized inverse was chosen. The gener- 
alized inverse gives a least squares solution to the overdetermined problem. The weighting 
was necessary because the Doppler sodar measurements are not as reliable as the surface- 
layer measurements. Weighting was accomplished by multiplying d and G by diagonal 
matrices, whose eigenvalues are interpreted as the variance associated with each datum. 
The surface-layer measurements were assigned a variance of one-tenth the Doppler sodar 
measurements. 

Since the modeling equations are nonlinear, the model parameters are determined by 
iterating the inverse problem. The model parameter estimate after the i + 1 iteration is 

m i+1 = (G)r 1 (d - di) + mi. (19) 


(17) 
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When there is good turbulent mixing (typically, the wind at 10 m > 3 m/s), the new 
method converges to a solution within ten iterations. However, the new method, like the 
older one, does not converge to a solution if the atmosphere is very stable. An example 
reconstruction using the new method, for 1548 on 8 March 1989, is shown in Fig. 3. The 
wind profile is on the left, and the temperature profile on the right. These profiles are 
typical of a convective boundary layer: the wind profile has a logarithmic shape, and the 
temperature decreases at a rate of one degree Celsius per 100 meters (the dry adiabatic 
lapse rate) outside of the surface layer. 

In Fig. 2, FFP predictions generated from the new method are compared with the 
older one. The sound speed profile is partitioned into layers at 5, 10, 20, 35, 50, 75, 100, 
125, 150, 175 and 200 m. The Doppler sodar data do not have a very significant effect 
on the predictions. In fact, the predicted levels are in even poorer agreement with the 
measurements than before. 

STOCHASTIC PROFILE RECONSTRUCTIONS 

Because the relationship between the meteorological profiles and the acoustic predic- 
tions is nonlinear, the •prediction from the mean sound speed profile is not necessarily the 
same as the mean prediction from the actual ensemble of profiles. Due to turbulent fluctu- 
ations, the actual sound speed profile varies about the mean. In this section an attempt is 
made to reconstruct an ensemble of realistic profiles; that is, an ensemble of profiles which 
include the fluctuating turbulent part. 

The reader may object to the neglect of the horizontal turbulent structure; whether a 
more realistic model for the horizontal structure would have much effect on the predictions 
remains to be determined. It should be kept in mind that the horizontal scale of turbulence 
is typically 100-500 m. (Ref. [3]) The vertical scale, however, is on the order of the height 
from the ground. For the purposes of acoustic predictions, it could be that a horizontally 
stratified atmosphere is a more realistic model than homogeneous turbulence. In any 
case, state-of-the-art modeling of homogeneous turbulence is discussed by Gilbert and 
Raspet [6], and the results of those authors are very similar to the results which will be 
given here. 

Stull [3], and Panofsky and Dutton [4], discuss parameterizations for the variances for 
atmospheric turbulence. The following axe in basic agreement with the equations presented 
in those two sources: For stable conditions, (<r^ = u'u')\ 

o u u» = 2.4, (20) 

<t t T, = 3.5. (21) 

For unstable conditions, 

/ 0 »i».\ 1 / 3 

Oulu , = (l2 - , (22) 
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/ lOzX -1 / 3 

or IT. = -2.0 (l - , (23) 

where r,- denotes the lower boundary of the stable layer, which caps the boundary layer. 

In the reconstructions of turbulence presented here, the fluctuations were ass um ed 
to have a jointly Gaussian distribution. The fluctuations at the various heights in the 
discretized profiles were assumed to have a correlation length given by 

t = kz. (24) 

It is a straightforward procedure to generate random numbers for a jointly Gaussian 
distribution on a computer, and the details of the technique will not be discussed here. 

As an example, the mean sound speed profile which was determined for 1548 on 8 
March 1989 is shown on the left of Fig. 4, and an ensemble of five stochastically-modified 
profiles is shown on the right. 

Predictions generated from the stochastically-modified profiles axe shown in Fig. 2. 
Agreement with the experimental data is improved, although the FFP still predicts a 
greater propagation loss than was actually measured. For each prediction appearing on the 
graph, ten profiles were generated. FFP predictions from each of these profiles were then 
averaged. Since the standard deviation of each ensemble of FFP predictions is about 10 dB, 
we can expect to be within about 3 dB (« 10 dB/\/IO) of the mean when ten profiles are 
averaged. Fluctuations resulting from this error are obvious in Fig. 2. Obviously, it would 
be desirable to average many more profiles, but this would increase the computational 
time unreasonably. For example, since the data shown on the figure required about 24 
hours of CPU time on a VAX Workstation, ten days of CPU time would be required to 
reduce the error to 1 dB. 


CONCLUDING REMARKS 

The modeling of meteorological profiles, and their use in predicting sound levels, was 
emphasized in this paper. Validity of Monin-Obukhov length similarity scaling was as- 
sumed in the modelling. The parameters needed for scaling were determined entirely from 
surface-layer measurements and remotely-sensed data, by performing a nonlinear inversion 
on mean meteorological data. The inversion technique presented in this paper, based on 
the generalized inverse, is flexible and works well when the problem is over- determined. 

Unfortunately, the inversion is not successful under the conditions of a very 
stable boundary. Monin-Obukhov model breaks down in this case, and no simple and 
accurate method for modeling profiles in the absence of good turbulent mixing have yet 
been developed. If realistic forward models can be developed for the stable boundary 
layer, the inversion technique presented in this paper could be easily modified, and should 
converge to a solution. 

The fast field program generates predictions for propagation in a horizontally- stratified 
medium. Therefore, only the effect of the vertical structure of turbulence on sound prop- 
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agation (and not the effect of the horizontal structure) can be modeled using the FFP. 
Nonetheless, since the horizontal scale of turbulence is typically much greater than the 
vertical scale, it still might be possible to model turbulence with sufficient realism using 
the FFP. In this paper, an initial attempt at turbulence modeling within the constraints 
of the FFP was made. Agreement between FFP predictions and acoustical measurements 
was improved, although at the cost of an order-of-magnitude increase in computational 
time. 
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LEVEL (DB RE= 0 AT 1 M) 


Rock Springs, 8 Mar. 1989 



Fig. 2 Comparison of experimental data and FFP predictions. The data (marked 

with a square) generally exceed the predictions. Predictions marked with a 
plus sign were generated from the sound speed profiles shown in Fig. 1. 
Predictions marked with a diamond include sodar wind profiles in the SSP 
reconstructions. Predictions marked with a triangle include sodar data and 
turbulence modeling. 
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Dote: 08 MAR 1989. Time: 1548 


Dote: 08 MAR 1989. Timet 1548 




wind speed (m/s) 


tempera tir'o (C) 


Fig. 3 Reconstruction of meteorological profiles for 1548 on 8 March 1989. Marked 

on the wind profile (left) are the sonic anemometer measurement at 10 m, 
and the sodar gates at 100, 150 and 200 m. Marked on the temperature 
profile (right) are the measurements at 1.9 and 8.9 m. 
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Date: 08 MAR 1989. Time: 1548 
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Fig. 4 On the right is the unperturbed sound speed profile for 1548 on 8 March 

1989. At left, an ensemble of five stochastically-reconstructed profiles, 
modeling the vertical structure of turbulence, is shown. 
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SUMMARY 

An FFP algorithm has been developed based on the method of Lee et al* for the 
prediction of sound pressure level from low frequency high intensity sources. In order 
to permit accurate predictions at distances greater than 2km, new correction 
procedures have had to be included in the algorithm. Certain functions, whose 
Hankel transforms can be determined analytically, are subtracted from the depth 
dependent Green's function. The distance response is then obtained as the sum of 
these transforms and the FFT of the residual k dependent function. One procedure, 
which permits the elimination of most complex exponentials, has allowed significant 
changes in the structure of the FFP algorithm, which has resulted in a substantial 
reduction in computation time. 

1 . INTRODUCTION 

Sound pressure levels at large distances from a point source close to the ground have 
been predicted using ray based procedures 1 in enhanced zones and residue 
calculations 2 in strong shadow zones. In the published literature (see references 1 and 
2) these predictions have been shown to be approximately valid. The errors for the 
predictions in the enhanced zone increase when ground reflections become important 
and when landing ray densities become small. In the shadow zone errors increase 
when the sound speed gradient becomes small. Both the above procedures are 
inaccurate or inoperable in the transition regions between shadow and enhancement. 

Since the publication of the first paper on the FFP 3 for atmospheric sound 
propagation this method has been increasingly used for sound pressure level 
prediction. This has largely occurred because the FFP can operate irrespective of 
whether there are shadow or enhanced or even mixed conditions present. Moreover 
the FFP can take proper account of ground reflection. 

The most widely known FFP algorithm, the CERL-FFP, stems from the method of 
Lee et al 4 , which was a development of the algorithm described in Raspet et al's first 
paper 3 . 

Starting from reference 4 we have reworked the analysis to enable us to produce our 
own FFP algorithm, structured in such a way that we could incorporate a variety of 
corrections in k space and thereby extend the range of validity of the result in the 
transformed (real) space. 


* Lee et al. J. Acoust. Soc. Am., 79, 1986, pp 628-634. 
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In section 2 of this paper a brief description is given of the development of our first 
prototype algorithm. A survey of previously published k spectral corrections applied to 
the first prototype is given in section 3. Section 4 describes our second algorithm 
which is based on a '7 averaging* procedure and section 5 describes a technique 
adopted for speeding up the computation. In the last section, 7, a comparison 
between the FFP and our other model's predictions is given for a realistic case. 

2. DEVELOPMENT OF THE FIRST FFP ALGORITHM 

2.1 Sign Convention 

Lee et al 4 replaced the system of atmospheric strata by an analogue electrical 
network. In our reworking of the model we found that this was not necessary and 
that retention of acoustic equations for pressure and particle velocity ensured greater 
clarity. In addition our analysis showed that great care had to be exercised with the 
signs used in the ladder calculation. 

Raspet 3 correctly drew attention to the need for different signs for the characteristic 
admittance dependent on the direction of the particle velocity, which gives the correct 
sign for the characteristic admittance of top stratum, Y co = + iy 0 /o)p 0 where y g is 
the propagation constant and p Q is the density and for the characteristic admittance of 
the ground Y c m = + 1/Z where Z is the usual ground impedance. 

None of the published papers on the FFP, including the most recent ones (see Franke 
et al 5 ), make it clear that the signs used in the ladder admittance calculation must be 
chosen in accordance with the direction in which the calculation is performed. The 
equation for calculating an admittance (Y new ) at the stratum interface nearer to the 
source from the admittance at the stratum's other interface (Yqj^) is 

Yrm W « — — ? nh + Y Q l d / Y c m (1) 

'new 1 cm \ 1 / 

1 + € (^old/^cm) tan h Y m C m 


where the characteristic admittance for layer m, is Y cm = e iy m /o)p m , 7 m being the 
propagation constant, p m the density and £ m the stratum thickness; c is +1 when 
working upwards and -1 when working downwards. We note that the e's cancel in 
(1) so that the same equation can be used whether working above or below the 
source interface. However in view of e in the top semi-infinite layer being set to +1 
the admittance calculated just above the source interface, Yfz^ - , has opposite sign to 
that calculated just below that height, Y(z s ) + . 

The z dependent Green's function at the source, Pfz^, is obtained from the known 
discontinuity in particle velocity at that location. Using our sign convention and 
noting the opposite signs of Y^)" 1- and Y(z s )“ a negative sign must appear in the 
denominator of P(z s ). 


P(z s > 


- (2i/a>p g ) 

Y(z s )+ - Y(z s ) “ 


( 2 ) 


Likewise performing the calculation of the Green's function at the detector, P(zq)* 
and retaining the above convention for t, the equation for calculation of p at the 
stratum interface nearest the source (P new ) from that at that stratum's other interface 
(P Q id) becomes 


Pnew “ P 


old | cosh 7m C m + £ sinh 7m C m [ J J < 3 > 


irrespective of whether the receiver is above or below the source. 


2.2 Ground Impedance 

The early FFP algorithms used the Delany Bazley 6 model for the ground impedance. 
Attenborough has suggested that his four parameter model 7 be used instead since this 
gives much smaller and more realistic impedance values at the low frequencies of 
interest in this study. 

2.3 Damping Coefficient 

The k spectrum calculated with the above ladder procedure has a large spike at the k 
value closest to co/c Q where c 0 is the sound speed for the top semi-infinite top layer 
and also a number of subsidary spikes at k values nearest to co/c m . These infinities 
produce errors when the Fourier transform is performed. A global damping is usually 
introduced, preferably by subtracting a small imaginary quantity ip from k, where p is 
typically 10” 4 . The transform is corrected for the effect of the damping by 
multiplying it by e/ 0 * 

This procedure effectively produces different damping effects on each spike dependent 
on their proximity to the nearest k sample point. The k sampling errors are 
therefore only partially removed. 

3. APPLICATION OF PUBLISHED k SPECTRAL CORRECTIONS 


3.1 Candel and Crance's k Shift Procedure 

Candel and Crance 8 proposed a method which ensured the k sample intervals were 
chosen so that the main peaks all occurred well away from the sample points. 
Thereby the transform errors arising from the presence of spikes in the k spectrum 
were substantially reduced. 

We set up an algorithm to readjust the k sampling interval, Ak, until the 5 main 
spikes were more than a limiting distance (Ak/10) from the nearest k sample point. 
The method requires a full ladder calculation for typically l/8th of the total k range 
centred on lqj = oj/c 0 for each setting of Ak. This can be time consuming when 
large numbers of strata are considered. 

The method worked well for large Ak values when the total number of k samples, N, 
was less than lk. For the larger N values required for ranges greater than 2km it 
became increasingly difficult to ensure the main spikes were distanced more than the 
above limiting value from the sample points. 

* The term k spectrum in this paper refers to F = k P(zjy) which is the function 
to be Fourier transformed . 



3.2 The Richards and Attenborough 9 Tail Remover 

The k spectrum very often has a non zero asymptotic value or tail as k goes to 
infinity. This tail invariably occurs when source and receiver are approximately the 
same height above ground. 

Introduction of a cut-off for the spectrum at a value k max produces large ripples in 
the transform which increase with range and with the proximity of k max to gj/c q . 3 
In this study we set k max at a value where the change in F was less than a 
prescribed limit, but this remained unsatisfactory. 

The k spectral tail can be removed by subtracting the function g(k) given by Richards 
and Attenborough 9 from F(k). 


g(k) - (A + e kz s) (1 - e ~ 6k ) 


(4) 


where Zg is the source height and A = \F (k max ) | and 5 is the derivative of 
\F (k)i at small k. g(k) has an exact Hankel transform which is added to the 
Fourier transform of F (k) - g (k). The above tail remover does not have the 
correct phase at small k and this produces a small error in the transform. 

Experiments with functions, with exact Hankel transforms which mimic the F (k) 
behaviour at small k and near the main spike are in progress. 

In Figure 1 the attenuation for a source and receiver 2m above an impedance ground 
at 50Hz is shown with and without the tail remover. 

4- MODIFIED FFP ALGORITHM USING 'y' AVERAGING 1 

For the large ranges of interest in this study very small Ak values must be used 
therefore the Candel and Crance method does not work. A novel technique for 
obtaining a k spectrum which can be more reliably transformed has been developed. 

The y averaging procedure is most easily understood by considering a single step in 
the admittance calculation for one stratum, index m, as described in section 2. 

Writing equation (1) in matrix form 


P 

U 


new 


P 


U 


old 


(5) 


where the matrix M m has elements which are functions of k and are usually evaluated 
at a single k sample value, say k r . 


M m- 


esinh y m Q m 
^cm 


cosh 7m C m 


( 6 ) 


where 7m -f k r 9 - u>Vc m 2 


I 


For Fourier integrals of smoothly varying functions F(k) (which have to be 
approximated by a discrete Fourier transform), the optimal sampling is at equidistant 
k r values. These F(k) samples can be assumed to approximate the contribution to the 
integral over the interval k r - Ak/2 to k r + /ik/2. In the presence of integrable 
infinities a different averaging process must be used. 

The simplest method follows from changing the integration variable in the above 
interval from k to 7; then elementary averaging leads to (6) with 7 m (k r ) replaced 

by 

Tm ” [ 7m O^r ~ + Tm Ohr + Ar/2) J/2 

We employ this averaging method for those intervals where 7 2 changes sign. 

Provided c m is real (no damping) the values of y m in the integrand are either pure 
imaginary or pure real and so are the 7 m values at all sample points except those 
where the above averaging is employed. Hence the hyperbolic functions in (6) can be 
replaced by real trigonometric or hyperbolic functions, the full complex functions being 
required at the above critical points only. A speed up by a factor of 3 was obtained 
by this procedure over that using (6) in its general complex form. This high speed 
algorithm cannot be used if artificial damping is present. 

5. k SPACE INTERPOLATION 

In order to obtain predictions out to large distances the number of k samples, N, for 
a given k,jj ax must be increased. As the major part of the computation lies in the 
determination of the z dependent kernel at each k value, we can achieve considerable 
improvements in calculation speed by interpolating F(k) where it varies smoothly. 

This applies to all regions of the spectrum lying outside the range spanned by oi/c m 
for all layers (typically k max /8). It was found that F(k) only needed to be evaluated 
every 8th point in the smooth region. 

For most cases linear interpolation of F(k) is adequate. However at very large ranges 
fluctuations occur in the transform due to the small discontinuities in the slope of the 
interpolated F(k). These errors can be reduced by using cubic interpolation. 

6. STRATUM QUANTISATION ERROR 

In Figure 2 the predicted attenuation above an impedance ground is shown for two 
different stratum configurations with the same small linear sound speed gradient (0.01 
s -1 ). For the solid curve all strata above 30m are taken as 27m thick and for the 
dotted curve as 54m thick. It is clear that the undulations in the dotted curve are 

much bigger than for the solid curve. This occurs because the FFP uses mid stratum 

sound speed averages for the whole of each stratum. This produces a stratum 

sampling error which gets worse the thicker the strata and the larger the sound speed 

gradient. There are two remedies: one is to use many very fine strata that increases 

the amount of computation, the other is to use linear sound speed variations in each 
stratum and use a modified procedure employing Airy functions 1 °- 

7. COMPARISON OF THE FFP PREDICTIONS WITH 

RAY/RESIDUE MODEL PREDICTIONS 

In Figure 3 the attenuation predicted by the FFP for a lapse-inversion-lapse 
meteorology is compared with that obtained by a hybrid method. This latter method 
combines the predictions from a ray based model for the enhanced regions with 
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those from a residue model for the shadow regions. The advantages of using the 
FFP are apparent in this figure. The residue model overpredicts the shadow 
attenuation and the ray model is restricted to giving predictions only in the region 
where there are ray landing points. Neither of these models can properly deal with 
the transitional region between shadow and enhancement. 
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ABSTRACT 


Weiner and Keast observed that in an upward-refracting atmosphere, the relative sound pres- 
sure level versus range follows a characteristic “step” function. The observed step function has 
recently been predicted qualitatively and quantitatively by including the effects of small-scale 
turbulence in a parabolic equation (PE) calculation. [Gilbert et al ., J. Acoust. Soc. Am. 87 , 
2428-2437 (1990)]. The present paper compares the PE results to single-scattering calcula- 
tions based on the distorted-wave Born approximation (DWBA). The purpose is to obtain a 
better understanding of the physical mechanisms that produce the step-function. The PE cal- 
culations and DWBA calculations are compared to each other and to the data of Weiner and 
Keast for upwind propagation (strong upward refraction) and crosswind propagation (weak up- 
ward refraction) at frequencies of 424 Hz and 848 Hz. The DWBA calculations, which include 
only single scattering from turbulence, agree with the PE calculations and with the data in all 
cases except for upwind propagation at 848 Hz. Consequently, it appears that in all cases ex- 
cept one, the observed step function can be understood in terms of single scattering from an 
upward-refracted “skywave” into the refractive shadow zone. For upwind propagation at 848 
Hz, the DWBA calculation gives levels in the shadow zone that are much below both the PE 
and the data. 


INTRODUCTION 

Weiner and Keast 1 and others 2,3 have observed that for sound propagation in an upward- 
refracting atmosphere, the relative sound-pressure level 4 versus range can be represented as a 
“step function” (Fig. 1). Recently the observed step function has been predicted qualitatively 
and quantitatively by parabolic equation (PE) calculations that include the effects of small- 
scale turbulence. 5 

Figure 2 shows gray-scale plots of the PE calculation without turbulence and with tur- 
bulence. The upward-refracted wave is called the “skywave.” 6 In the plots with turbulence 
the skywave is still present although it has been noticeably modified by turbulence. For a 

* This work was supported by the Office of Naval Research 
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source and receiver near the ground, the average relative sound pressure level inside the sky- 
wave (region 1 in Fig. 1) is approximately 0 dB (spherical spreading) with or without turbu- 
lence. However, the region below the skywave (region 3 in Fig. 1) is dramatically affected by 
turbulence. Without turbulence a deep shadow zone is predicted by the PE calculation. With 
turbulence, sound is scattered from the skywave into the shadow zone, producing a relative 
sound pressure level that is fairly uniform on the average. The region between the skywave 
and shadow zone (region 2 in Fig. 1) is a transition region. The horizontal extent of region 
2 is a strong function of the strength of upward refraction. It is evident that, for a gray-scale 
plot with turbulence in Fig. 2, a horizontal “cut” through the plot at a particular receiver 
height will give a step function. 

Although the gray-scale plots of the PE calculations give a good qualitative picture for un- 
derstanding the step function, the PE calculations do not allow a simple physical explanation 
of the observed quantitative behavior of the relative sound pressure level versus range. For ex- 
ample, what causes the observed constant relative sound pressure level (spherical spreading) 
at long ranges? In the present paper a simpler calculation is presented which is based on sin- 
gle scattering out of the upward- refracted skywave. The simpler calculation, which uses the 
distorted- wave Born approximation (DWBA), 7 is compared to the PE calculation and to the 
data of Weiner and Keast. The objective is to gain insight into the physical mechanisms that 
produce the observed step function. 


I. THEORY 
A. Atmospheric model 

We want to compare the DWBA calculations to the PE calculations reported in Ref. 5. 
Hence we use the same atmospheric model as in Ref. 5 and assume that the effects of tur- 
bulence can be adequately represented by small fluctuations in the index of refraction. The 
total index of refraction is thus written as a steady deterministic part n d (R) plus a fluctuating 
stochastic part fi(R,t) where R — (x,y,z) and t is time. With this approximation for turbu- 
lence, the wavenumber is given by 


k(R,t) = k 0 [n d (R) + fi(R,t)] , (1) 

where k 0 is a reference wavenumber, n d = 1, and fi • Cl. The deterministic part of the index 
of refraction n d is assumed to vary only with the height above the ground z. It was computed 
from a logarithmic sound-speed profile of the form 


( Co + a ln(z/d), z > z 0 
l Co + a ln(z 0 /d) z < z 0 


( 2 ) 


where Co = 340 m/s, z 0 = 0.01 m, and d = 6 x 10 -3 m. The refraction parameter a is —.5 
m/s for crosswind propagation (weak upward refraction) and -2.0 m/s for upwind propagation 
(strong upward refraction). The deterministic parameters were chosen to fit the short-range 
data of Weiner and Keast. 1 
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In the calculations reported here, the fluctuation fi is assumed at any instant of time to be 
a function of R = (x,y, z). In Ref. 5 the turbulence model was two-dimensional so that sound 
propagated only in the x — z plane. As will be shown, for computing the average sound pres- 
sure level using a single scattering approximation, the two atmospheric models are equivalent. 

In Ref. 5 the stochastic wavenumber in Eq. (1) was used directly to calculate “snapshots” 
of the acoustic field. Here, we want to compute average levels so we need the autocorrelation 
function for \i. The autocorrelation function is defined by 


C(S) = {fi(R + S)fi(R)) , (3) 

where {) denotes an ensemble average over many realizations of //. (We assume an ensemble 
average and time average are equivalent.) For small-scale turbulence near the ground, C(5) 
can be approximated by a Gaussian autocorrelation function of the form 

C(S ) = ^2 e -(S|//|+52//|+5?//2) 5 (4) 

where fi 0 is the root-mean-square value of //, and l x , l y , and / 2 are the correlation lengths in the 
x, y and 2 directions, respectively. In numerical calculations isotropic turbulence was assumed 
(L = ly = lz = 0 - The input values for fi 0 and l (fi 0 = 1.42 x 10 -3 , and / = 1.1 m) were taken 
from measurements reported by Daigle. 8 


B. Ground impedance model 

The ground was modeled as a flat, locally reacting plane with an angle- independent com- 
plex impedance. Impedance values were obtained from the empirical formulas of Delaney and 
Bazley 9 using an effective flow resistivity of 300 rayls/cm. The resulting impedance values 
were 7.19 + *8.20 and 5.17 + *5.57 at 424 Hz and 848 Hz, respectively. 


C. Distorted-wave Born approximation (DWBA) calculations 

We consider a point source with angular frequency u in a turbulent atmosphere. At a par- 
ticular instant in time the solution for a point source (Green’s function) in a turbulent atmo- 
sphere satisfies 


V 2 G'(R, R') + k 2 0 {n d + y) 2 G{R , R') = -AwS(R - fr) , (5) 


where R! is the source location, and R is the receiver location. In the absence of turbulence 
(fi = 0) the Green’s function Go is given by 


V 2 G 0 (R, #) + kln 2 d G 0 (R,R!) = -4t tS(R - Rf) 


(6) 
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where Go is a refracted wave (i.e., a distorted wave) if n<* varies with height. In this section we 
shall consider both an undistorted spherical wave and a wave distorted by upward refraction. 

An integral equation for G which goes to Gq in the absence of turbulence can be written as 


G(R,R') = G 0 (R,R') + i- ! Go{R,R!')Sk\R'')G{k',R')d 3 R" , (7) 

47T J 


where Sk 2 = &o( n d + l 1 ) 2 ~ ^o n d — 2/ifc, o , since = 1 and /< < 1, 

Equation (7) allows us to write the total solution, G, as the solution in the absence of tur- 
bulence, Go, plus an integral which gives the turbulent contribution. However, since the un- 
known G appears in the integral, Eq.(7) is as difficult to solve exactly as is the original differ- 
ential equation. When Sk 2 is “small enough” the full solution G that appears in the integral 
can be approximated by G 0 . The approximation G = Go is generally known as the “Born 
approximation”. When Go is a refracted wave the approximation is often called the “distorted- 
wave Born approximation” or “DWBA.” 7 

Writing the turbulent contribution as SG , using G = Go, and Sk 2 = 2/j.k q we have 

SG=^J G 0 (R, R")n(R")Go{R", R')d 3 R" . (8) 

We want to calculate the time average of |G| 2 , which we assume is the same as an ensemble 
average. Denoting the average value of \G\ 2 as (|G| 2 ) and assuming a random phase between G 
and SG we have 

<l<?l 2 > = {| Go \ 2 ) + <|«?| 2 ) , (9) 


where 


d & I 2 ) = / G'oA R'")GlA", #)<„(#'>(«")} (10) 

X Go(R, R")Go(R", R')d 3 R , "d 3 R" . 

Now (fi(R'")n(R")) = C(S) where C is the autocorrelation function and S = R ,n — R” . Trans- 
forming from the variables (R", R'") to the variables (R" , S) gives 

X Go(R, R")Go{R", R')d 3 s d 3 R" . 
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For the sake of illustration we first consider an undistorted spherical wave in free space, i.e., 
we take Gq(R,R") = exp(ifc&|./2 — R"\)I\R — /2"|. For |/?'| = 0 (source at the origin) and 
\R"\ » |S|, the Green’s function Gq(R" + S,R' = 0) can be approximated as 


G„(R" + S, K = 0) = 


(12) 


where the k = koh, and n = (R" /\R"\). Similarly the Green’s function Go(R, R" + S) can be 
approximated as 


.ik 0 \R-R"\ 

Go(R,R" + S) * =; — e 

| R-R”\ 


(13) 


where k' = k 0 h ', and h 1 = (R — R")/\(R — i?")|. With these approximations for the free-space 
Green’s functions, we have 


(I | 2 ) = ^ [ -4 ,-U— s C0)<PS<PR" 

Vl 1 ' 4 tt 2 J I R" \ 2 \R- R" |2 v ' 


We now define a scattering function a(q) as 


a(q) = J e i? S C(S)d 3 S , 


where q — k' — k. Then Eq. (14) becomes 


(14) 


(15) 


(| 8G | 2 ) = P- [ -J— <r(g) ■ 1 „ d*R" 

' 4tt 2 J I R" p kh, \R- R" 12 


(16) 


In Eq. (16) we can identify l/|.ft"| 2 as the sound intensity I tnc incident on the scattering vol- 
ume and 1/| R — R" | 2 as the scattered intensity I scat that reaches the receiver. Hence we can 
write 


(ISGf) = ■^JlUR")a(q)I, M (^')<PR" . (17) 

Equation (17) has a useful physical interpretation (see Fig. 4). The average intensity of the 
sound reaching the receiver from a particular volume of space is proportional to the product of 
the incident intensity reaching the volume, the scattering strength of the volume, and the scat- 
tered intensity. The Appendix gives an analytic result for Eq. (17) for small-angle scattering. 
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In order to take upward refraction into account we use the following prescription: In Eg. 
(17) we replace the incident intensity 7, nc = 1/|7?"| 2 and the scattered intensity 7, nc = 1/|7? — 
R" | 2 with Ii nc = \Gpe{R") | 2 and I scat = \Gpe(R,R")\ 2 , respectively, where Gpe denotes the 
Green’s function without turbulence computed using the parabolic equation method described 
in Refs. 5 and 10. Writing the integrals in Cartesian coordinates we have 


(I 6G | 2 ) = J | Gpe(R'') 2 | G pe (R,R ") | 2 

x exp [ i(q x S x + q y S y + q z S z )] C(S X , S y , S z )dx" dy” dz" dS x dS y dS z 


(18) 


Since Gpe is azimuthally symmetric and I72"! >> |5|, we neglect the y-dependence in Gpe and 
integrate exp (iq y S y ) over y and obtain the 6-function result obtained in the Appendix. (See 
Eqs. (A2) - (A4).) Also we set q x = 0, as in the Appendix, and obtain 


(I 1 2 > = J H / *"(* - *") I Gpe(R") PI G pe (R, R") p (ig 

xe iq * s >C{S x ,Q,S z )dS x dS z dx"dz" . 

Using the general Gaussian autocorrelation function in Eq.(4) for the integrations over S x and 
S z , we have 


(I SG P) = I X " {X - x ") I a PE (x'\ z") p (20) 

x | Gpe(x, Z] x", z") I 2 e- q * l *t 4 dx"dz" 

In the parabolic equation method the quantity actually solved for is ^(r, 2 ), where Gpp(r, z) — 
[exp(ifcor) / y/r ]^(r, z), and r = yjx 2 + y 2 . Since the integral is now two-dimensional, we can 
set y to zero and let r = x. Then in terms of ^(x, z) we have 


(I SG | 2 ) 


y 0 k 0 l x l z f T / // , /x ,2 1 lTf / 12 

— Y x — J I 2 ) II - x i z -< z ) I 

xe-W A dx"dz" . 


( 21 ) 


In the numerical calculations we assumed the turbulence to be isotropic with a correlation 
length /. Hence we have finally 


(I SG | 2 ) 


^ J | x, ** p| *(*-*»;*,*") P 

xe-&' 2 l 4 dx"dz" , 


where rr ;/ goes from zero to x, and z n goes from zero to infinity. 


(22) 
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II. NUMERICAL CALCULATIONS 

A. Comparison of DWBA calculations with PE calculations and with experiment 

In Fig. 3 the DWBA calculations and the PE calculations are compared to each other 
and to the data of Weiner and Keast. The data is for octave bands of random noise between 
300-600 Hz and 600-1200 Hz, respectively. In both the DWBA and the PE calculations, the 
frequency was taken to be where fi and / 2 are the lowest and highest frequencies, re- 

spectively, in the octave bands considered. Section I gives the parameters for the atmospheric 
model and the ground impedance model used in the calculations. Note that the DWBA calcu- 
lations and the data are for the average relative sound pressure level while the PE calculations 
are a “snapshot” of the relative sound- pressure level and not the average level. However, the 
trend in a particular PE calculation is generally fairly close to the average level predicted by 
the corresponding DWBA calculation. 

The DWBA calculations, which include only single scattering from turbulence, give a good 
approximation to the average PE levels in all cases except for upwind propagation at 848 Hz. 
For upwind propagation at 848 Hz, the DWBA prediction deep in the shadow zone is much 
below both the PE and the data. 


B. Discussion of numerical results 

Some of the features of the curves in Fig. 3 can be understood in a straightforward 
way using the DWBA calculations. The deterministic (no turbulence) part of the Green’s 
function is G 0 and the stochastic part due to turbulence is SG. Near the source (regions 1 and 
2 in Fig. 1) we have | Go | 2 ^ (| SG | 2 ) while at long ranges (region 3 in Fig.l) we have just the 
reverse. Consequently, near the source, the shape of a given curve for relative sound pressure 
level versus range is governed by the deterministic sound-speed profile so the level is essentially 
what one would obtain from a calculation without turbulence. 

Since we have | Go | 2< C (| SG | 2 ) at long range, the relative sound pressure level is due al- 
most entirely to scattering from turbulence. In order to understand the long-range behavior of 
the curves in Fig. 3 we must make a more detailed analysis than was required at short range. 
As shown in the Appendix, the contribution to the relative sound pressure level from turbu- 
lence scattering in free space (with no refraction) diverges as the logarithm of the range. We 
expect similar behavior even with upward refraction over a finite impedance plane. Consider 
the situation in Fig. 4 where we have a scattering volume with an incident intensity 7, nc and a 
scattered intensity I aca t- The sound intensity incident on the scattering volume is proportional 
to 1/r 2 where r is the horizontal range to the receiver. The scattering volume itself is propor- 
tional to r 3 . The scattered intensity reaching the receiver from the scattering volume, like the 
incident intensity, is proportional to 1/r 2 . For a fixed scattering angle, the average scattered 
intensity from the volume is thus proportional to (1/r 2 ) x (r 3 ) x (1/r 2 ) = 1/r. Hence, as 
shown in Eq. (A9), we expect the relative sound pressure level to increase as the logarithm 
of the horizontal range. This behavior at long range is seen in the DWBA calculation for 
crosswind propagation (weak upward refraction) at 424 Hz. When there is significant upward 
refraction the height of the scattering volume is not proportional to the range but increases 
more rapidly than linearly with range. As a result, the scattering angle is not fixed but 
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increases with increasing range. Since the scattered intensity is reduced at larger scattering 
angles, the relative sound pressure level versus range is “flattened” so that a nearly constant 
relative sound pressure level is reached at long ranges. Because of the flattening effect caused 
by an increasing scattering angle, a nearly constant relative sound pressure level is seen in the 
DWBA calculation for upwind propagation (strong upward refraction) at 424 Hz. A similar 
behavior is seen for crosswind propagation (weak upward refraction) at 848 Hz. The flattening 
effect with weak upward refraction at 848 Hz is apparently due to the greater sensitivity to the 
scattering angle at the higher frequency. 

The DWBA calculation for upwind propagation (strong upward refraction) at 848 Hz falls 
off in the shadow zone much more rapidly than the PE calculation and the data. The major 
computational difference between the two calculations is that the PE calculation includes mul- 
tiple scattering while the DWBA calculation does not. Hence the disagreement indicates that 
for upwind propagation at 848 Hz, multiple scattering is important. This interpretation is sup- 
ported by a gray-scale plot for this case which shows the skywave greatly modified by turbu- 
lence so that the approximation G = Go is not valid. 


III. SUMMARY AND CONCLUSIONS 

We have compared distorted-wave Born approximation (DWBA) calculations to parabolic 
equation (PE) calculations and to the data of Wiener and Keast. In all cases except one, the 
DWBA calculations, which include only single scattering, predicted the step-function behavior 
of the relative sound pressure level versus range seen in both the data and the PE calculations. 
The important conclusion to be reached is that, in the presence of upward refraction, single 
scattering can give a relative sound pressure level that does not diverge as the logarithm of the 
range but rather is nearly constant at long range. Hence, in all cases except one, the observed 
step function can be understood in terms of single scattering from an upward- refracted sky- 
wave. 

For upwind propagation (strong upward refraction) at 848 Hz, the DWBA calculation 
grossly underestimated both the data and the PE calculation. In this case, the single scat- 
tering approximation G = Go was not valid in the skywave. To accurately predict multiple 
scattering of sound into the shadow zone, one must have a good predictive model for sound 
propagation in the skywave itself. Hence, it would be valuable to have measurements not only 
for the sound scattered into the refractive shadow, but also for the sound field in the skywave. 
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APPENDIX: ANALYTIC EXPRESSION FOR FREE-SPACE SCATTERING FROM 

ANISOTROPIC TURBULENCE 

To see the general behavior of Eq. (17) we can consider weak small- angle scattering in free 
space . For weak small-angle scattering we can use the Born approximation and obtain an an- 
alytic result for anisotropic turbulence. 

For small angle scattering we can let 1/|R"| 2 = 1/x" 2 and 1/|R — R " | 2 = 1 /(x — x") 2 . We 
could integrate Eq. (17) directly using a particular autocorrelation function such as a Gaus- 
sian. However, to obtain a more general result that does not assume any particular autocorre- 
lation function, it is convenient to return to the form in Eq. (14) which is written in terms of 
the autocorrelation function C(S X , S y , S z ), 


<l iG l 2 > -£/ 


1 

2 


C(S X , Sy , S z ) 


p i(<lxS x -\-qySy+qzSz) 

(x - x") 2 

dx" dy" dz" dS x dSy dSz . 


For small angles we can approximate q x as 


(Al) 


1- = >* (i + -TTZ-y 


X — x" 


=k o 


[x"(x — x") 


Similarly, 


q y = k 0 


x 


[x"(x — X W )J 


We now consider the integral over z": 


(A2) 


(A3) 


[ + °° e iSzkolx/x"(x-x")]z" dz » _ 2ir S 
J — oo 


ko- 


x H (x - x") 
2n x"{x — x") 
ko x 


S(S Z ) . 


(A4) 


The integral over y" gives a similar result with S z replaced by S y . Inserting the results 
from integrating over z" and y" into Eq. (Al) and integrating over S z and we have 

(| SG | 2 ) = § / e’ , * 5l C’(S' x ,0, 0) dx" dS x . (A5) 

We could integrate over S x and define a special scattering function 
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(Too = f C(S X , 0, 0) exp (iq x S x )dS x . However, we are considering small-angle propagation. Hence 
q is almost perpendicular to the propagation direction and we can therefore set q x to zero. 
Thus, integrating in the region between the source and receiver we have 

(I SG | 2 ) = 


where 

Coo 

For anisotropic turbulence having a Gaussian autocorrelation function (See Eq.(4)) we ob- 
tain 


U 2 too tx 

/ C(S X , O, O) dS x / dx" 
X Z J-oo Jo 


k 2 

K o n 
— » 

x 


(A6) 


[°° C{S X , O, 0)dS x . 

J — oo 


(A 7) 


(| SG | 2 > = ^ 4 kl l./x . (A8) 

Thus for weak small-angle scattering in free space, the scattering due to turbulence falls off 
inversely with the range and depends on the correlation length only in the direction of propa- 
gation. Note that, written in terms of the relative sound pressure level (RSPL), the contribu- 
tion from turbulence scattering diverges as the logarithm of the range. 

RSPL = 101og w (x 2 (l SGI 2 )) (A9) 

= iOlogjo (V^o lx) + 101og 10 (x) . 
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RELATIVE SOUND 
PRESSURE LEVEL (dB) 



Fig. 1. Characteristic “step” function for the relative sound pressure level versus range for 
sound propagation in an upward-refracting atmosphere. (From Gilbert et al. 5 ) 
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Fig. 2. Gray-scale plots of relative sound pressure level versus height and horizontal range 
for a non-turbulent atmosphere and a turbulent atmosphere. The frequency is 424 Hz, 
and the source height is 3.7 m (12 ft). Parameters for the atmospheric model and ground 
impedance model are given in the text. (From Gilbert et al. 5 ) 
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Fig. 3. Relative sound pressure level versus horizontal range for a refracting, turbulent at- 
mosphere. The connected circles are the data of Weiner and Keast. 1 The solid lines are 
parabolic equation results from Ref. 5. The connected x’s are distorted-wave Born ap- 
proximation calculations. Parameters for the atmospheric model and ground impedance 
model are given in the text. 






Scattering 

Volume 




Fig. 4. Schematic representation of scattering from turbulence. The quantity /, nc is the 
average intensity incident on a particular scattering volume, and I scat is the average 
scattered intensity. The total scattered intensity is obtained by integrating over the vol- 
ume between the source and receiver. 
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WAVE PROPAGATION THROUGH RANDOM MEDIA: 

A Local Method of Small Perturbations based on the Helmholtz Equation 

Ralf Grofie, Universitat Oldenburg, Postfach 2503, D-2900 Oldenburg, FRG 


1 Introduction 

Propagation of sound through the turbulent atmosphere is a statistical problem. The randomness of the 
refractive index field causes sound pressure fluctuations. Although no general theory to predict sound 
pressure statistics from given refractive index statistics exists, there are several approximate solutions to 
the problem. The most common approximation is the parabolic equation method. Results obtained by 
this method are restricted to small refractive index fluctuations and to small wave lengths. While the first 
condition is generally met in the atmosphere, it is desirable to overcome the second. This paper presents 
a generalization of the parabolic equation method with respect to the small wave length restriction. 

2 Parabolic Equation Method 

For the small wave length limit 1 the Helmholtz equation can be converted into a parabolic form (main 
propagation direction e z ) f\/: 

( 2 ^' + S? + & + *'■«) 0M = O M 

k = wave number; p = refractive index deviation; ip = complex sound pressure 

The refractive index deviation p is considered as a random function. Therefore equation (l) is a stochastic 
differential equation and the sound pressure tp becomes a random function, too. Stochastically this equation 
is nonlinear, e.g. it contains a product of two random variables. For this reason the stochastic parabolic 
equation cannot be solved exactly. Further approximations are necessary. Several mathematical tools 
were applied to remove the stochastical nonlinearity, i.e. path integrals /2/, functional derivatives /3/, 
perturbation expansion methods /!/. Despite approximations used in the calculations looking different, 
the results are all the same. 

The physical meaning of all these approximations becomes evident in the local method of small perturbations 
/ 1/ . By this method the scattering volume is divided into slabs perpendicular to the main propagation 
direction. Each slab is chosen as thin as required by the validity limit of the first order perturbation 
expansion term (single scattering approximation, Born approximation). This distance clearly depends on 
the strength of the refractive index fluctuations. If these fluctuations are sufficiently small, the slabs are 
much thicker than one correlation length of the random medium. Therefore the slabs can be regarded as 
uncorrelated. Based on both assumptions - small refractive index fluctuations and a small wave length 
compared to the correlation length - the statistical independence of subsequent slabs can be proofed 
mathematically. Wave propagation through random media is described here as a Markov process. 

*The wave length must be small compared to the size of a typical inhomogeneity of the medium. In statistical terms this 
size is expressed by the correlation length of the refractive index autocorrelation- function. 
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As a consequence of the Markov property slabs of finite thickness are no longer necessary. This results in 
a differential equation for the mean sound pressure which is linear in the stochastical sense as well 2 /!/: 

< t/>(r) >= 0 (2) 

( 3 ) 


< V'(r) >= ^o(r) exp{-az} (4) 

tpo — free propagated incident wave 

The mean complex sound pressure decreases exponentially while the wave propagates through the random 
medium. This is an effect of decorrelation of the sound wave due to phase fluctuations. Different members 
of the statistical assembly interfere destructively because of their different phases. 

The validity of result (4) - and of any other result obtained by the parabolic equation method - is restricted 
first by the validity of the parabolic wave equation (1) and second by the validity of the Markov assumption. 
Necessary conditions are the smallness of refractive index fluctuations and the smallness of the wave length 
(compared with the correlation length). In the next section the first condition is also assumed to be 
true. The small wave length assumption, however, will be replaced by the weaker condition of negligible 
backscattering. This will lead to a generalized Markov process and, consequently, to generalized results 
with respect to the wave length - correlation length ratio. 

3 Generalized Local Method of Small Perturbations 

While generalizing the parabolic equation method the main idea of the local method of small perturbations 
will be retained. The refractive index fluctuations are assumed to be small enough to justify the application 
of a single scattering approximation within a distance A z in the scattering volume which is large compared 
to the correlation length. Again the scattering volume is divided into slabs of this size. Therefore subsequent 
slabs are uncorrelated as well. But contrary to the parabolic equation method the contribution of one slab 
will not be calculated from the parabolic equation but from the Helmholtz equation: 

(a + k 2 ( 1 + M (f))) tp(r) = 0 (5) 

A = Laplace operator 

Neglecting backscattering 3 yields a difference equation for the mean sound pressure: 

< ip(p n ,nAz) >= (g 0 + £) < (n - l)Az) > (6) 

2 It is possible to derive similar equations for the higher statistical moments of V* by the parabolic equation method, too. 
Only the first moment equation and its solution are presented here in order to compare them with the generalized results 
derived in section 3 of this paper. The decorrelation coefficient a is calculated for an exponential autocorrelation function for 
comparison, too. 

3 Neglecting of backscattering is not an assumption of the parabolic equation method but a consequence of the small wave 
length limit. The parabolic equation implies that there are small scattering angles only and no backscattering at all. 


/ d d 2 d 2 \ 

( al + + v + j 


a = 


<n 2 > k 2 i 


l = correlation length 
Equation (2) can be solved easily: 
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p — (x, y) ; n = number of slab 

A A 

Go is the integral operator of the homogeneous Helmholtz equation (p = 0, free propagation) and S is 
an integral operator for the scattering within one slab. Since double scattering is the lowest order non- 

A 

vanishing term, the kernel of S contains the autocorrelation function of the refractive index field: 

S {r nt r) < tp(f) >= -k 4 Jd*? Jd?r G(f n , f')G(f' ,r) < > < H*) > ( 7 ) 

G = Greens function of the Helmholtz equation 

The solving of equation (6) is somewhat different from that of the related parabolic equation problem. 
Here no wave length approximation helps to calculate the integrals. But for the case of a homogeneous 
refractive index autocorrelation function, equation (7) becomes a convolution product. It can be Fourier- 

transformed with respect to the variable p, and this operation turns the operators S and Go into simple 

ni rs/ 

functions S and Go. In the Fourier representation equation (6) reads: 

<V> (k, nAz) >— + s^j <ip ( K , (n - l)Az) > (8) 

k = 2-dimensional spatial frequency 

By Fourier-transformation the ^-integrations are already performed and only the z-integrations are left. 
They can be performed, too, if the z-dependence of the medium autocorrelation function is known. For 
the sake of simplicity this dependence is assumed to be exponential 4 . After z-integration the scattering 
contribution of one slab is seen to be proportional to the slabs thickness Az. Therefore equation (8) can 
be written as (5 = s Az): 


<V» (/c,nAz) >= + 8 A <«/> (/c, (n — l)Az) > 

The effect of all slabs is obtained by iteration: 

<^> (k, nAz) >= ^G 0 + s A z^ (k,0) 

Regarding a sufficiently large number of slabs yields (z = nAz): 

»v /v 

<t/> ( k,z ) >=V> 0 («.*) exp{a (/c)z} 

For the special case of an isotropic exponential autocorrelation function s becomes: 


( 9 ) 


( 10 ) 


( 11 ) 


< ft 2 > k 4 
4 oaj(aj - a) 


( 12 ) 


o = Vk 2 — /c J , aj = -^(k + i/l) 2 - k 2 (13) 

To compare (11) and (12) with the corresponding parabolic equation method results (4) and (3), the incident 
wave V»o is assumed to be a plane one. Then the Fourier-transformation of (11) results in: 

4 For more complicated functions the s-integration can be performed numerically. 
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(14) 


< i K*) >= V>o(0 ex p{ * (°) z ) 


5(0) = 


< m 2 > * 2 / (i« - J fe 2 / 2 ) 

4 (1 + fc 2 / 2 ) 


(15) 


The real part of 5 (0) describes the decorrelation of the sound wave. It is a more general expression than 
(3) - only in the small wave length limit they are equal. Because of s(0) being a complex number, a second 
effect is predicted by this method, which cannot be seen in the parabolic results. The imaginary part 
is a stochastic correction to the wave number k due to an increase of the mean propagation distance in 
the random medium. Only in the small wave length limit the scattering angles are small and the mean 
propagation distance corresponds to the z-extension of the scattering volume. 


4 Conclusions 

A generalized form of the local method of small perturbations has been presented in this paper. Working 
directly from the Helmholtz equation instead of the parabolic equation the small angle scattering method 
was replaced by a forward scattering method . By this method only one result wets derived here: The first 
statistical moment of an incident plane wave scattered by a very weakly statistical homogeneous random 
medium with an exponential autocorrelation function. This result shows corrections to the corresponding 
parabolic equation method result. 

It is possible to apply the method to more complicated problems, i.e. a difference equation for the second 
statistical moment can be derived by the same idea. 

If the medium fluctuations become stronger, the thickness of one slab decreases. The slabs might be thicker 
than the correlation length, but not as much as assumed before. Then correlations between two slabs have 
to be taken into account. This leads to difference equations which connect statistical moments not only 
from one slab to the following, but to the next following, too. 

The most valuable advantage of this method might be its suitability for numerical calculations. For any gi- 
ven medium autocorrelation function the scattering function s can be obtained by Fourier-transformation. 
The incident wave is also Fourier-transformed. Then the algorithm given by equation (9) is applied iterati- 
vely until the desired propagation distance is covered. The final result is obtained by Fourier-transformation 
again. 

The author wishes to thank Prof. K. Haubold, Prof. V. Mellert, Dr. M. Schultz-von Glahn and A. Sill for 
the valuable discussions during the whole work. 
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Abstract 

A system has been designed to provide an assessment of noise levels 
that result from testing activities at Aberdeen Proving Ground, MD. The 
system receives meteorological data from surface stations and an upper air 
sounding system. The data from these systems are sent to a meteorological 
model, which provides forecasting conditions for up to three hours from the 
test time. The meteorological data are then used as input into an acoustic 
ray trace model which projects sound level contours onto a two dimensional 
display of the surrounding area. This information is sent to the 
meteorological office for verification, as well as the range control office, and 
the environmental office. To evaluate the noise level predictions, a series of 
microphones are located off the reservation to receive the sound and 
transmit this information back to the central display unit. The computer 
models are modular allowing for a variety of models to be utilized and 
tested to achieve the best agreement with data. This technique of prediction 
and model validation will be used to improve the noise assessment system. 

Introduction 

The U.S. Army has an active testing program for munitions and 
weapons located at Aberdeen Proving Ground, MD (APG). The results of 
these tests can cause high sound levels to impact on the local community. 
This problem has existed for a long period of time, but recently it has 
become more acute because of the development of the local communities 
adjacent to the Proving Ground. APG is actively engaged in a number of 
different programs to alleviate the noise problem. One of the approaches to 
mitigate noise complaints is to be able to better indicate when conditions 
could enhance the sound propagation at long distances due to the 
atmospheric structure. As a result of these concerns, the Noise Assessment 
and Prediction System (NAPS) was proposed utilizing sensors, models, and 
computers to predict the noise levels that might be encountered at an 
off- range site as a result of a particular test. 
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Objective and Approach 


The objective of the NAPS development is to create an operational 
system for predicting noise intensities based upon present and forecasted 
diurnal meteorological conditions. The reason for specifying diurnal 
conditions is to limit the meteorological model to only those conditions that 
change with solar input. The meteorological model will not be used to 
forecast synoptic conditions, passage of fronts, etc; synoptic scale conditions 
accounted for utilizing standard weather forecast techniques and tools. 

The NAPS development approach is to modify and adapt existing 
acoustic and meteorological prediction models for this noise prediction 
problem. The aim is to be able to have these models operate on a PC 
located at the meteorological office, providing the information to the 
various range offices responsible for testing. In order to provide timely 
information to the users, the results of the noise predictions will be made 
available to users every 15 to 20 minutes. Users will have an assessment of 
the current conditions and how they may vary within the next three hours. 

These criteria required that both the acoustic and meteorological 
models be compact with short execution times in order to meet the 
required specifications. Therefore, it was initially decided to utilize a 
standard ray trace model 1,2 with modifications for its use to make 
predictions at APG. In the same vein, a 1-D planetary boundary layer 
model 3,4 was chosen and incorporated into NAPS. 

In the development of NAPS, it was decided to utilize a SODAR which 
provides wind averages and the occurrence of wind shears at 15 minute 
intervals. The SODAR measurement coupled with other meteorological 
measurements from an instrumented mast at the test site and upper air 
data from a radiosonde would provide the required input data for the 
predictive meteorological model. To better aid the meteorologist and range 
personnel in determining the propagation conditions, the data are 
assimilated from the different sensors, processed through the various models 
to provide displays of the meteorological profile, the ray trajectories, or the 
contour of sound intensities overlaid on a terrain map of APG and 
presented at each users office. These computer displays aid in making the 
test scheduling and GO/NO GO decisions. 

The next step in the development process is the evaluation of the 
system to determine its performance and fine tune the system to an 
operational performance level for use on a daily basis in support of the test 
programs at APG. 

Data Requirements 

The data requirements and operations for NAPS consist of a 
radiosonde measurement at 8:00 AM, to provide information on the 
atmospheric conditions up to 5 km. The number of radiosonde flights 


depends on the synoptic conditions, ranging from a minimum of one release 
for no synoptic changes during the day to a number determined by the 
meteorologist monitoring the changing synoptic conditions. The radiosonde 
provides vertical profiles of temperatures, winds, and relative humidity from 
near surface to 5 km. The vertical profile can detect for occurrences of 
temperature inversion and wind shear conditions which can cause the sound 
to be refracted to the surface rather than being refracted upward. 

Wind conditions within the planetary boundary layer (PBL) (surface to 
1-2 km), whose height varies diurnally, are monitored by a SODAR. The 
SODAR is a remote sensor which provides 15 minute averages of winds and 
wind shears to approximately 600 meters. This permits a continual update 
of the atmosphere near the surface; the part of the PBL subject to the 
greatest changes during the progression of the day. As mentioned 
previously, there are two-meter meteorological masts located at each of the 
test locations. These measure surface temperature, winds, humidity, 
pressure, and solar radiation. In the future, plans call for adding a 
ten-meter mast. This would permit measurements of meteorological 
parameters at the two-meter and ten-meter levels. The two and ten-meter 
configuration will enable meteorologists to utilize similarity theory and other 
techniques to model the surface layer meteorological conditions. Again, the 
vertical extent of the surface layer varies and is dependent upon the solar 
radiation input, the type of surface and wind speed. 

Data from the various sensors will be continually monitored by the 
meteorologist to ensure the accuracy of the observations. The data is then 
entered into NAPS to provide an assessment of the present conditions and 
how these conditions vary under the influence of diurnal and terrain 
conditions. Once the meteorological data is verified, it is provided as input 
into the acoustic propagation model (ray trace) to calculate ray trajectories 
and noise intensity contours. These are again examined by the 
meteorologist to verify that the predicted intensities at the different 
locations are reasonable and agree with the meteorological conditions. The 
meteorologist, after verifying the data is consistent, can now release the 
data to the range personnel to assist them in making a decision about 
upcoming tests. 

Acoustic Models 

NAPS provides an estimate of peak noise intensity from a blast source 
at ground level in all directions from the blast source. An essential feature 
of the model is its ability to account for variations in meteorological 
conditions in the calculation of sound propagation. In performing noise 
intensity estimates, acoustic ray traces are generated each 5° (or multiples 
of 5°) in azimuth, over a sufficient range of elevation angles to define the 
focusing and shadow regions in the area surrounding the blast. The NAPS 
model accounts for spherical spreading, absorption 5 , focusing, shadow 
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zones, reflection of rays from water, interference of multiple rays arriving at 
the same location, the directional asymmetry of a blast, and the terrain 
elevation. Essential model inputs are the vertical temperature, vector wind 
speed, humidity structure of the atmosphere, the blast charge weight, blast 
location, and blast height. 

Meteorological Model 

The acoustic model requires input from meteorological sensors or 
meteorological parameters derived from a meteorological prediction model. 
NAPS is required to be able to predict favorable propagation conditions 
and when conditions are not favorable for the test. To accomplish this, a 
1-D planetary boundary layer model was acquired and adapted to operate 
on a PC. The initial meteorological conditions are input to the model 
utilizing the surface meteorological data at the firing site, SODAR data 
from the adjacent location and the upper air data from the standard 
morning sounding or a sounding closest in time to the test period. The 1-D 
model generates a vertical profile of temperature, vector wind speed, and 
humidity from the surface to the top of the boundary layer. 

The measured data from the sensors are merged into a single wind and 
temperature profile at the site. The profiles with additional meteorological 
measurements and the geostrophic winds at 850 mb are used as input into 
the 1-D Planetary Boundary Layer model. The model is then used to 
predict meteorological profiles at one hour intervals up to three hours in 
time. These profiles are then used as input into the ray trace model to 
predict acoustic intensity levels resulting from a particular test and firing. 

System Description and Operation 

The various components and sensors comprising the NAPS system are 
shown in figure 1. Data from the various sensors are linked to the PC in the 
Meteorological Station by either hardline or RF link. The data is ingested 
into the PC, evaluated, and then sound contours are predicted for a 
particular test. The meteorological data, both measured and predicted from 
the model, are transferred to the Range Control Office, where it is used as 
input into the same acoustic ray trace model as being run at the 
meteorological station. This permits the Range Control Office to share the 
same information that is available at the Meteorological Station. A view of 
the data flow in the system is shown in figure 2, where the data are used as 
input into the meteorological model. Examples of this output are shown in 
figures 3-6 which are the wind speed, direction, temperature and speed of 
sound, respectively. From this point, the data are input into the acoustic 
model with output from the acoustic model shown in figures 7 and 8 as ray 
trajectories generated at given azimuths. Also displayed are the speed of 
sound profiles showing the atmospheric structure which causes the rays to 


be refracted either upward or downward. Sound level contours are 
generated by utilizing the ray traces at 5° increments from 0° to 360°. In 
addition, for post analysis and system evaluation, data from the off-post 
microphones are collected and put into the computer for inclusion in the 
graphic display for archiving with the meteorological data. 

To demonstrate how NAPS would operate, meteorological data are used 
as an input into the acoustic model which produces the sound level contours 
shown in figure 9. These contours are generated from measured data and 
indicate what the sound intensity levels would be at the present time. The 
contours are elongated and could result in some loud noises impacting upon 
the local community. The next step is to determine how the situation might 
change in the next several hours. Prediction of the meteorological 
conditions for one and three hours later are made by the meteorological 
model and inserted into the acoustic model with the results shown in 
figures 10 and 11. In figure 10, one hour later, the changes in the contours 
are appreciable, with the overall contour shape becoming rounder. Three 
hours later, there are some changes in the contours, but these are not as 
significant compared to those showing the change from present time to one 
hour later. Reviewing the data, as it becomes available, indicates that the 
test might be delayed to an hour until conditions for testing have improved. 

System Evaluation 


The situation at APG is excellent for evaluating meteorological and 
acoustic models since the sound source characteristics 6,7 and locations are 
known; and there are a large number of atmospheric sensors located 
throughout APG. To verify the complaints and provide comparisons for 
NAPS, a noise monitoring system is used to provide measurements of the 
propagated sound levels. Figure 12 is a map which shows the location of 
the meteorological and acoustical sensors on and off APG. The asterisks in 
figure 12 indicate the seventeen microphone sites which are located off 
APG. These are set up to operate at a threshold of 108 db. When the noise 
exceeds this level, it is recorded and transmitted with the time of 
occurrence to a computer at range control and from here it is transmitted 
to the meteorological station. There are five surface meteorological masts 
sited on APG and three that are located off APG to the east, west, and 
north of the Proving Ground. In addition, there are two SODARS located 
approximately 12 miles apart which provide winds and wind shear heights; 
these are shown by the open circle. Upper air soundings are made at the 
meteorological station which is also indicated by an open circle located 
adjacent to the SODAR at the north end of APG. These sensors then 
provide detailed data on the meteorological conditions at APG, and the 
microphone monitoring system provides sound level intensities from those 
tests that exceed a level of 108 db. 
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This configuration of sensors and sources can provide a system for 
evaluating the acoustic predictions made by the ray trace model as well as 
the predictions of meteorological conditions made by the meteorological 
model. 

As mentioned earlier, preliminary evaluations have been made of the 
NAPS prediction capability. An example of this is shown in figure 13 where 
there is a fair agreement between the predicted sound level contours and 
those levels measured by the microphones. 

It is planned to evaluate the performance of the NAPS over a minimum 
of an annual cycle, since there are seasonal periods when the occurrence of 
high intensity off range are greatest. To be able to capture the required 
data, a NAPS data base management system is being developed. Figure 14 
is a diagram of this system. There are two major parts to the system; one 
is located at APG, the NAPS operational site, and the other site is located 
at WSMR* which is the prototype development site. The WSMR site will 
be used to test and evaluate the software and hardware before integration 
into the operational NAPS at APG. The data base will consist of data 
obtained at both sites, which have markedly different environments from 
each other. In the case of WSMR, the environment is a desert one, with 
low humidity, higher temperatures and greater solar radiation. The APG 
site is more a continental maritime site situated on the Chesapeake Bay. This 
site would be more humid, with lower temperatures and less solar radiation due 
to the presence of clouds, vegetation, and inclement weather. It will be 
interesting to compare similar type data from each of the sites. By 
analyzing data from both sites, it may be possible to gain further insight to 
local variations at each of the sites, thereby making the utilization of NAPS 
at other locations easier. 

Summary 

The NAPS was developed to predict sound level intensities resulting 
from testing at APG. NAPS utilizes standard in-situ meteorological sensors 
in addition to remote sensors. A ray trace acoustic model and a 1-D 
planetary boundary layer model are used to predict sound propagation 
conditions out to three hours based on the meteorological model. A data 
base is being developed to capture the acoustic and meteorological data 
and to utilize this data on evaluating and improving the sound intensity 
predictions. The data will include data from at least a years period to 
insure that NAPS has been evaluated and utilized under a variety of 
diurnal and seasonal conditions. After a thorough evaluation, the NAPS 
will become an operational system. The information learned by putting this 
type of operation at APG can then be used in installing the NAPS at other 
sites that may be having a noise problem which could be mitigated by 
taking into account the effects of the atmosphere on acoustic propagation. 
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FIGURE 5. COMPUTER DISPLAY OF TEMPERATURE PROFILE. 
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FIGURE 6. COMPUTER DISPLAY OF SOUND SPEED PROFILE. 
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FIGURE 7. RAY TRACE WITH SOUND SPEED PROFILE. 
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FIGURE 8. RAY TRACE WITH SOUND SPEED PROFILE. 
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DB CONTOUR PLOT 



FIGURE 9. SOUND LEVEL CONTOURS DERIVED FROM METEOROLOGICAL MEASUREMENTS. 



FIGURE 10. SOUND LEVEL CONTOURS DERIVED FROM PREDICTED METEOROLOGICAL 


CONDITIONS ONE HOUR LATER. 




FIGURE 11. SOUND LEVEL CONTOURS DERIVED FROM PREDICTED METEOROLOGICAL 
CONDITIONS THREE HOURS LATER. 



FIGURE 12. ACOUSTIC AND METEOROLOGICAL SENSOR LOCATIONS AT ABERDEEN 
PROVING GROUND. 
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13. PREDICTED SOUND LEVEL CONTOURS COMPARED TO MEASURED SOUND 
LEVELS FROM MICROPHONE SYSTEM. 
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FIGURE 14. NOISE ASSESSMENT AND PREDICTION SYSTEM DATA BASE CONFIGURATION. 
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ABSTRACT 

According to ray theory, regions exist in an upward refracting atmosphere where no sound should be 
present. Experiments show, however, that appreciable sound levels penetrate these so-called shadow 
zones. Two mechanisms contribute to sound in the shadow zone: diffraction and turbulent scattering of 
sound. Diffractive effects can be pronounced at lower frequencies but are small at high frequencies. In 
the short wavelength limit, then, scattering due to turbulence should be the predominant mechanism 
involved in producing the sound levels measured in shadow zones. No existing analytical method 
includes turbulence effects in the prediction of sound pressure levels in upward refractive shadow zones. 

In order to obtain quantitative average sound pressure level predictions, a numerical simulation of the 
effect of atmospheric turbulence on sound propagation is performed. The simulation is based on scattering 
from randomly distributed scattering centers ("turbules"). Sound pressure levels are computed for many 
realizations of a turbulent atmosphere. Predictions from the numerical simulation are compared with 
existing theories and experimental data. 




INTRODUCTION 

Solar heating of the ground produces strong temperature gradients in the air just above the surface of 
the Earth. Since the speed of sound is proportional to the square root of the temperature, sound will 
follow upwardly curved paths in every direction from a source. The stronger the temperature gradients 
involved, the shorter the distance to what is properly called a shadow zone, since no direct or reflected 
rays can penetrate into this region. Figure 1 is an illustration of an upward refractive shadow zone where 
the source is at a height h s and the radius of curvature of the limiting ray is The edge of the shadow 
zone is delineated by the so-called limiting ray which grazes the ground. 

In a similar fashion, sound traveling upwind is curved upwards due to the strong wind gradients 
near the ground and a shadow zone is also formed. In the case of wind, however, the effect is not 
isotropic because of the vector nature of the wind velocity and the rays are actually bent downward for the 
sound propagating downwind. 

Two mechanisms contribute to the magnitude of the sound levels measured in shadow zones: 
diffraction and the turbulent scattering of sound. Pierce 1 describes the solution for a linear sound speed 
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gradient in terms of residue series for the pressure in the shadow zone. He examines the cases of a hard 
boundary and a pressure-release surface and gives approximate solutions when the source and receiver are 
above the creeping wave layer height, defined as (Rg/2 where R^. is the radius of curvature of the 

limiting ray and k 0 is the wavenumber. 

Daigle et al. 2 made use of the above two approximate solutions in an effort to fit the data they 
collected over an asphalt airport runway and over a grass-covered strip near the runway, the latter 
approximating a pressure-release surface at frequencies greater than 500 Hz and the former approximating 
a hard boundary. They found that the hard boundary data was well explained by Pierce's approximate 
solution for that case and that the data up to 1000 Hz over the grass-covered ground was satisfactorily 
explained by Pierce's approximate solution for a pressure-release surface. 

The approximate solution leads to large errors in the effective source levels for sources close to a 
pressure-release or finite impedance ground, as was the case with Daigle's data. A complete discussion of 
this problem can be found in the paper by Raspet and Franke. 3 In a later paper. Berry and Daigle 4 used 
the complete residue series solution and again compared the above data. They found that the data at 250 
Hz still agreed well with the predictions of diffraction theory. But the data was well under-predicted by 
diffraction theory at 500 Hz and especially at 1000 Hz. The predictions from the full residue series 
solution are shown in Fig. 2, which is a reproduction of Fig. 13b of Ref. 4. 

The role played by atmospheric turbulence in the insonification of shadow zones has escaped 
analytical formulation. In an effort to obtain a quantitative estimate of the extent to which atmospheric 
turbulence raises the sound levels in a shadow zone, Gilbert et al. 3 used a parabolic equation method to 
numerically simulate sound propagation in a turbulent atmosphere. They compared their predictions for 
upward refracting conditions with experimental results of Wiener and Keast.6 The numerical predictions 
involved the calculation of the sound pressure magnitude for a particular realization or "snapshot" of 
turbulence, while the results of Wiener and Keast were expressed in terms of average sound pressure 
levels. Nevertheless, Gilbert et al. were able to duplicate the apparent range independence of excess 
attenuation characteristic of the experimental data at ranges as great as 1 km. 

In this paper, we present the average sound pressure levels in an upward refractive shadow zone 
predicted by a scattering center based numerical simulation. The main features of the numerical solution 
are reviewed and the modifications necessary to adapt it to an upward refractive atmosphere are discussed. 
Sound levels are computed for over 500 realizations of the turbulent atmosphere. Predictions from the 
numerical simulation are then compared with experimental data taken by Daigle et al. 2 

MAIN FEATURES OF THE NUMERICAL SIMULATION 

Although the details of the numerical simulation were given in an earlier paper* the main features are 
repeated here so that the reader may have a better idea of the type of calculations involved. Following the 
model of de Wolf, 7 we construct an ensemble of isotropic, irrotational scattering centers which we call 
"turbules." If |i is defined as the change from unity of the index of refraction, a given turbule is assigned 
the refractive profile 

* Walton E. McBride, Henry E. Bass, Richard Raspet, Kenneth E. Gilbert, "Scattering of sound by 
atmospheric turbulence," submitted to J. Acoust. Soc. Am., Feb. 1990. 
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( 1 ) 


[L(r,s) = qie' r/S 

where < 7 j is the value of |1 at the center of the spherically symmetric turbule and s is the 1/e contour of 
the scattering center and can be considered to be its effective size. The value of < 7 , and the probability 
distribution of turbule sizes depend in general on the particular functional form chosen for the correlation 
function of the fluctuations of the index of refraction. If the correlation function is chosen to have the 
Gaussian form, 


) = ^ e ' r/L > (2) 

where <p 2 > is the variance and L is the correlation length, then the size spectrum is a delta function 
implying that all the turbules have the same size. 



The value of < 7 ,- for this particular form of the correlation function is given by 


(3) 



(4) 


and is inversely proportional to the turbule number density p^. An upper value of p^j of about half the 
overlap density is necessary so that the turbules will be separate entities. With single scattering, sound 
scattered from a particular turbule reaches the receiver downfield with negligible scattering by other 
turbules located between that particular turbule and the receiver. From Eq. (4) the product < 7 ,- 2 p N is a 
constant whose value depends on the independently measured micrometeorological variables <p 2 > and L. 
There is, therefore, a certain latitude in the value of p^. Decreasing p N will result in a greater value for 
l< 7 ,l. Although a lesser number of turbules result from a decrease of p N , the predictions of the numerical 
simulation are statistically steady as the turbule number density is decreased from an upper limit of half the 
overlap density. 

Initially in the development of the simulation, the first Bom approximation to scattering was used to 
determine the scattering effect of each turbule. In practice, the evaluation of the scattering integral is 
performed by assuming that both source and receiver are far away from the scattering region; thus first 
order terms in the phase are sufficient. Because some of the turbules in the numerical simulation are close 
to the source or receiver, second order terms in the phase were kept in the scattering integral. The total 
pressure at a receiver downfield is, then, the sum of the direct and scattered spherical waves. For only one 
turbule in free space, this is: 
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(5) 


p(R) = 


R 



j^2 s 3 e «'*(*> n-o,) 


1 ^ e -CJfcV/4 


where 


C = (l-cos ©o) 2 + sin 2 ©o 
and 



( 6 ) 


(7) 


In the above equation, k is the wavenumber, s is the effective size of the turbule, © 0 is the angle 
between the incident and scattered directions, R is the distance between source and receiver, while r st is the 
distance between source and turbule center, and r^ is the distance between turbule center and receiver. 

Note that the usual Bom scattering term is recovered when a = 0 and the first term in Eq. (6) is 
dropped. Even with this improved evaluation of the Bom scattering integral, the distance from turbule to 
source or receiver cannot be less than about twice the radius of the turbule. Consequently, "buffers" of a 
turbule's diameter were placed in front of the source and receiver where no turbules were allowed. 

The numerical simulation using the first Bom approximation to scattering was then compared to 
theoretical expressions due to Karavainikov^ for the log-amplitude and phase variances of the pressure 
fluctuations. It was found that good agreement was reached whenever the wave parameter D (=R/&L 2 ) 
was greater than 1. As shown in Fig. 3, the log-amplitude variances as predicted by Karavainikov are 
independent of frequency when D < 1, a result also obtained by Bergmann using geometrical optics. 

In an effort to reach better agreement in the geometrical optics region, the Rytov approximation used 
by Karavainikov was incorporated into the numerical simulation. The Rytov method consists of 
approximating the field at the receiver by 


p R (f) = p 0 (r)eV R ®; 


( 8 ) 


whereas the Bom approximation is written: 


p B (r) = p 0 (f) + ? B (?X 

There is a simple relationship between the first Rytov and first Bom approximations: 


(9) 


Po(?) 


(10) 


and it was, therefore, a simple matter to incorporate the first Rytov approximation into the numerical 
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simulation. The results are shown in Fig. 4 and good agreement is obtained throughout the range of the 
wave parameter D. 

As can be seen from the above comparisons with Karavainikov's analytic curves, the first Rytov 
approximation is superior to the first Bom approximation for an unbounded medium without refraction. 
When refractive conditions are introduced, however, the formation of shadow zones becomes possible. In 
the shadow zones, p 0 is 0 and the first Rytov approximation cannot be used. Recourse must be made to 
the first Bom approximation and the wave parameter D must be greater than 1 for the numerical simulation 
to be valid in accordance with the results of Fig. 3. 

The next step is the inclusion of the ground. An immediate consequence of the existence of a 
boundary is the presence of three additional paths by which sound can propagate to the receiver. There 
now exist four single scatter paths that connect the source and receiver: 

1 . source-turbule-receiver, 

2. source-turbule-ground-receiver, 

3 . source-ground-turbule-receiver, 

4 . source-ground-turbule-ground-receiver. 

The last three paths all interact with the ground and, therefore, a model of the effect of the ground on the 
sound wave was also included in the numerical simulation. 

The algorithm proceeds as follows. Values of <(i 2 > and L are given from independent micro- 
meteorological measurements. From these, the value of s and are obtained using the above equations. 

A scattering space, which will enclose thousands of turbules, is defined with buffers in front of the source 
and receivers of widths equal to about the diameter of a turbule. The turbules are assigned positive or 
negative signs for their value of q v The sound pressure at the receiver is calculated for this particular 
arrangement of turbules, and the result is referred to as a realization. Then each turbule is given random, 
small increments in its Cartesian coordinates. The sound pressure at the receiver is recalculated, resulting 
in another realization. The process is repeated for as many realizations as are necessary for the statistics to 
stabilize. We have found that 500 realizations are sufficient. Average sound pressure levels can then be 
obtained from the 500 stored values of the sound pressure. It should be mentioned that any other desired 
statistical quantity can be obtained, such as structure and correlation functions, as well as the variances of 
the log-amplitude and phase fluctuations. 

The inclusion of a sound speed gradient requires two modifications: the rays are now curved and the 
value of the wavenumber k is no longer constant along a ray. In order to obtain a closed form solution 
for the equation describing the rays, a linear sound speed gradient was assumed. As is well known, a 
consequence of this assumption is that the ray paths are arcs of circles. 

Because of the curvature of the ray paths, each path must be tested to see whether the source and 
each turbule can be joined together, as well as each turbule and the receiver. If either segment of the total 
path cannot be linked, that particular turbule's contribution is discarded. It was found that about 15% of 
the turbules were eliminated in this way for the experimental data to be described later. 

The last correction necessary is the calculation of an effective wavenumber k e for each ray path. 

This required the computation of the length of each path, as well as the travel time along that path. Their 
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ratio gave an effective sound speed c e along that path, and the effective wavenumber was then given by 
oVc e . 

COMPARISON TO DATA 

In order to implement the numerical simulation, the statistical properties <|i 2 > and L of the 
turbulent atmosphere and the impedance of the ground are required. The former was given in the article by 
Daigle et al. 2 as <p 2 > = 6 x 10' 6 and L = 1.6 m. The impedance had to be approximated because the 
article mentioned above did not specify a particular impedance model. To estimate the impedance, a 
residue series solution developed by Raspet and Franke^ was used to match the curves of Fig. 2 as closely 
as possible at all three frequencies. The particular impedance model used was a four parameter model 
developed by Attenborough.^ A shape factor ri of .750, a shape factor ratio Sj - of .875, a porosity Q of 

.675, and a flow resistivity a of 330 cgs rayls give the results shown in Fig. 5. 

The numerical simulation was performed with the above parameter values for 500 realizations. Rms 
sound pressure values were computed and divided by the pressure doubling plus attenuation factor 
(2 e'QK/R) as was done in Daigle's presentation of his experimental data. The source was given a height 
of 0 m, and six receiver positions were used 230 m downrange at heights of 0.25 m, 0.50 m, 1.0 m, 2.0 
m, 4.0 m, and 7.0 m in order to sample the vertical behavior of the sound pressure levels in the shadow 
zone. Because we are forced to use the first Bom approximation when the receiver is in a shadow zone 
(pQ = 0), the possible values of the wave parameter D must be checked to see that they will be greater than 
1. Since D = R/fcL 2 with R = 230 m and L = 1 .6 m, it is sufficient to check the value of D for the greatest 
frequency of interest. With 1000 Hz, k ~ 18.4 and, therefore, D = 5 which is well above the minimum 
value of 1. The sound pressure levels predicted by the numerical simulation (expressed in dB) are 
compared to Daigle's experimental data in Fig. 6 for the three frequencies involved. Notice the similarity 
in shape which the curves representing the turbulent scattering contribution share with the curves that are 
typical of diffraction theory predictions. As can be seen by comparing Figs. 2 and 6, it appears that for 
this experiment the contributions from turbulent scattering and diffractive effects are about equal at 250 
Hz. However, turbulent scattering becomes the major contributor at 500 Hz. At 1000 Hz, turbulent 
scattering is the predominant mechanism behind the increased sound pressure levels measured in the 
shadow zone. 


CONCLUSION 

Quantitative predictions for the average sound pressure levels in a refractive shadow zone have been 
presented. Good agreement was reached with experimental data collected by Daigle et al 2 in a shadow 
zone caused by temperature gradients. It seems that the use of a simple linear sound speed gradient is a 
good approximation to the real sound speed profile directly above the ground for the moderate ranges 
involved in this study. At the longer ranges investigated by Gilbert et al. it was necessary to use a 
logarithmic sound velocity profile to obtain accurate predictions.^ The relative contributions of diffraction 
and turbulent scattering have been examined and graphically displayed. The dominant mechanism which 
dictates sound levels in shadow zones at higher frequencies is scattering due to turbulence. 
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Fig. 1 Shadow zone formation for a sound speed profile that decreases linearly with height. 



Fig. 2 Comparison of measured sound levels (points) with predictions based upon diffraction into the 
shadow zone, taken from Ref. 4. Solid circles are for 1000 Hz, triangles for 500 Hz, and 
diamonds for 250 Hz. 
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Fig. 3 Comparison of the numerical simulation using the first Bom approximation (O for phase variances; 

X for log-amplitude variances) with Karavainikov's analytical results ( for phase variances; 

for log-amplitude variances). 



Fig. 4 Comparison of the numerical simulation using the first Rytov approximation (o for phase 

variances; X for log-amplitude variances) with Karavainikov's analytical results ( for phase 

variances; for log-amplitude variances). 
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Fig. 5 Comparison of data with the prediction of diffraction theory using the residue series solution. Data 
points are described for Fig. 2. 



Relative dB Level 


Fig. 6 Comparison of the prediction from the numerical simulation with experimental data taken in 
shadow zone. Experimental data is the same as in Fig. 5. 
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Abstract 

Acoustic measurements made in the atmosphere have shown significant 
fluctuations in amplitude and phase resulting from the interaction with 
time varying meteorological conditions. The observed variations appear to 
have short term and long term (1-5 minutes) variations at least in the 
phase of the acoustic signal. One possible way to account for this long term 
variation is the use of a large scale wind driven turbulence model. From a 
Fourier analysis of the phase variations, the outer scales for the large scale 
turbulence is 200 meters and greater, which corresponds to turbulence in 
the energy-containing subrange. The large scale turbulence is assumed to 
be elongated longitudinal vortex pairs roughly aligned with the mean wind. 

Due to the size of the vortex pair compared to the scale of our experiment, 
the effect of the vortex pair on the acoustic field can be modeled as the 
sound speed of the atmosphere varying with time. The model provides 
results with the same trends and variations in phase observed 
experimentally. 
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Effects of Large Scale Wind Driven Turbulence 
on Sound Propagation 


Introduction 


Random fluctuations in the acoustical index of refraction in the 
atmosphere is the result of the presence of turbulence. These random 
fluctuations in the acoustical index of refraction results in fluctuations of 
the amplitude and phase of an acoustic wave. The variations in the 
amplitude and phase show changes occurring over two different time 
scales. The short term variations correspond to turbule sizes on the order of 
1 meter, while the long term variations seem to correspond to turbule sizes 
on the order of 100 m and greater. 

The aim of this work was to develop a descriptive model for large scale 
wind driven turbulence and the effects of large scale turbulence on the 
sound field. The model will describe the shape and horizontal and vertical 
wind velocity profiles for the turbulence. Due to the size of the turbules in 
relation to the experiment conducted, a simple phase model was developed 
to perform phase variation calculations using the results from the large 
scale turbulence model. 

Atmospheric Effects 

Before the model for the large scale wind-driven turbulence is presented, 
lets first examine the dynamics of the atmosphere. In discussing the details 
of air flow, it is convenient to consider the atmosphere to be divided into a 
number of horizontal layers (figure 1). The region in which the atmosphere 
experiences surface effects through vertical exchanges of momentum, heat, 
and moisture is called the planetary boundary layer (PBL) or is sometimes 
referred to as the friction layer. Panofsky and Dutton 1 defines the depth of 
the PBL, h, as the thickness of the turbulent region next to the ground 
which is also called the mixing layer. Another height used to describe the 
thickness of the PBL in the daytime is the height z, of the lowest inversion. 
Actually, h tends to be roughly 10% larger than z; because the lowest part 
of the inversion is still turbulent, partly because of overshooting from 
below, partly because there is often strong wind shear in the inversion. 

The lowest part of the PBL is called the surface layer. In this layer, the 
characteristics of turbulence and the vertical distribution of mean variables 
are relatively simple. There is no precise definition of the surface layer. 
Qualitatively, the surface variations of vertical fluxes can be ignored. 
Typically, the fluxes are large at the surface and decrease to zero near the 
top of the PBL. 
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The main problem is calculating the height of the lowest inversion z,-. 
This value is important since it represents the largest size an inhomogeneity 
can be in the atmosphere. According to Panofsky and Dutton 1 , the 
horizontal wind speed fluctuations are related to z,- by 


u. 


(12 - 0 . 5 - p -) 1/3 

■bmn 


( 1 ) 


where u. is the friction velocity and L mo is the Monin-Obukhov length. If 
variations in the horizontal wind speed are due to purely mechanical 
turbulence, an alternate formula for it. can be used for z > z 0 : 


uk 

ln(z/z 0 ) 


( 2 ) 


where k is the von Karmon constant (0.4), u is the horizontal wind speed at 
height z, and z Q is the roughness length. Substituting equation (2) into 
equation (1) and solving for z,- results in 

Zi = 2Z mo [12 - (^) 3 /n 3 (z/z 0 ) ] (3) 


This provides the height of the lowest inversion in terms of Monin-Obukhov 
length, the fluctuation of the horizontal wind speed, and the roughness 
length. The Monin-Obukhov length can be estimated using tables 1 and 2 
knowing the surface wind, incoming solar radiation, and the roughness 
length (for table 2, the roughness length was 0.05 meters for the 
experiments 1 ). 


Experimental Procedure and Data Analysis 

A series of line-of-sight propagation measurements were made over 
relatively flat open farm land. A run consisted of an eight minute record of 
signals received simultaneously at five transverse microphones mounted one 
meter above the ground and one microphone mounted near the source for a 
reference (figure 2). The sound source was driven by a tape with a 
prerecorded signal consisting of a mixture of eight tones centered at one 
octave spacings beginning at 62.5 Hz. This geometry is similar to the 
geometry Daigle 2,3 used in his experiments. 

The meteorological data was collected using a series of three-cup 
anemometers and temperature probes at four heights; 3, 10, 30, and 110 ft. 
The data acquisition system provided a five minute period of wind speed, 
wind direction, and temperature as well as the maximum and minimum 
values during the five minute period. Measurements of the fluctuating wind 
speed and temperature data were also made using the techniques outlined 
by Johnson 4 . 

The Fourier transform of the amplitude and phase variations contains 
the spectrum of the fluctuations of the sound field due to turbules present 
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in the atmosphere. The spectral peaks are related to the scale of turbulence 
L by “Taylor’s hypothesis of frozen turbulence” which relates the temporal 
and spatial turbulence scales by 5 

L = ut (4) 

where u is the mean wind speed and r is the characteristic time associated 
with the temporal measurements. Taylor’s equation can be rewritten as 

L = H (5) 

v 

where v = 1 /r. Calculations of L show the different scales of turbulence 
present in the atmosphere during the experiment. Figure 3 is for a run 
where the wind speed is low. The spectrum shows several peaks which 
represent the different scales of turbulence present in the atmosphere for 
that run. Figure 4 is for a run where the wind speed is high. The only 
spectral peak present is one at a low frequency. This implies that the only 
scale of turbulence which is affecting the phase is on the order of a few 
hundred meters in size. 

Some caution must be noted here about this type of analysis. The 
location of the low frequency peak may be a result of insufficient frequency 
resolution due to the length of the sample analyzed. A longer time sample 
might shift the low frequency peak to even lower frequencies. 

The Fourier transform for the amplitude variations were also computed. 
There is not a spectral peak for the amplitude at the low frequency end of 
the spectrum. Large scale variations in the atmosphere cause changes in 
the sound field resulting in refractive variations instead of a scattering 
process as in small scale turbulence. 


Large Scale Turbulence Model 


The first problem is to obtain, from experimental measurements, a clear 
idea of the structure and motion of the turbulence. From now on, frequent 
references will be made to ‘eddies’ of the turbulent motion, a word intended 
to describe flow patterns with spatially limited distributions of vorticity 
and comparatively simple forms. Since the experimental data consists of 
point measurements, the identification of eddy types must be by informed 
guesswork followed by measurements designed to confirm the guess. 

According to Tennekes, 6 there appears to exist in all turbulent shear 
flows more or less distinct large eddies with relatively long lifetimes. 
Townsend was the first to investigate the structure and dynamics of these 
large scale vortices. 7 Townsend was struck by the fact that in all turbulent 
shear flows he knew, the eddy viscosity K m , nondimensionalized by 
appropriate length and velocity scales, turned out to be a number that is 
relatively independent of the flow considered. Townsend hypothesized that 
the large eddies must be responsible for this universality. According to 
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Townsend, the eddies are elongated longitudinal vortex pairs in the 
boundary layer, roughly aligned with the mean flow, figure 5. 

The lifetimes of the eddies are greater than the length of time for data 
runs discussed here. For this analysis, they will be considered to be 
“permanent” . Note that secondary circulations cause local regions of 
horizontal convergence near the surface. Those regions are the sites of 
vigorous turbulence production rates, and may be responsible for the 
generation of most of the Reynold stress in the boundary layers. Tennekes 
concludes that the eddies are capable of relatively long lifetimes because 
the mean shear is an adequate source of energy. 

If a stream function /(x, z) for a particular arrangement of eddies is 
known, there are several parameters of the eddy system which can be 
calculated. Stream functions are a type of function which describe the 
streamlines in a flow. Streamlines are regions where the velocity vectors of 
the fluid are tangent at a particular instant. The velocity distribution of 
the eddy can be calculated using 8 


and 


u(x,z) 


dfjx, z) 
dz 


(6) 


v(x,z) 


dfjx, z) 

dx 


( 7 ) 


where tt(x, z) and v(x, z) are the horizontal and vertical wind speeds 
respectively. The functional form of the stream function which represents 
an eddy pair is 7 


f(x,z) = A[cos(lx) + e ,2 7“*]e «“ 2r2 (8) 

where A is a constant specifying the intensity of the eddy pair, 
a 2 r 2 = a 2 x 2 + ct 2 z 2 , l is the characteristic wavenumber of the eddy pair, 
and ck x and a z are the horizontal and vertical wavenumbers for the eddy 
pair. The coordinates (x,z) are relative to the center of the eddy pair. 
Townsend uses a characteristic wavenumber for the eddy pair of na x . 

Using equations (6) and (7), the horizontal and vertical wind speed are 

u(x, z) = Ba 2 z z[cos(lx) + e -,2 /° r2 ]e - « or2r2 (9) 

and 

u(x, z) = — B{2lsin(lx) + a 2 x[cos(/x ) + e -,2 7 a2 ]}e“« 0(2r2 (10) 

where B = A/ 2. Figure 6 is the horizontal wind speed versus height for x 
= 0 m, a x = 0.0043 m -1 , a z = 0.0087 m -1 , and B = 2000 m 2 /s. The 
negative height refers to a vertical position below the center of the eddy 
pair. Figure 7 is the horizontal wind speed versus range for z = -150 m 
using the same parameters as in the previous figure. The negative range refers 
to a horizontal position to the left of the center of the eddy pair. 
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For this work, the size and intensity of the eddy pairs were determined 
from meteorological data taken in the field. The standard deviation of the 
wind speed was calculated using a TSI hot wire anemometer. Using the 
standard deviation of the wind speed, the roughness length, the wind speed 
at height z, and estimating the Monin-Obukhov length from table 2, the 
height of the lowest inversion layer z,- is calculated using equation (3). This 
provides a maximum height of the eddy pair. The Fourier transform of the 
phase variation provides an estimate of the lower limit for the horizontal 
extent of the eddy pair. For the data analyzed, the average of the wind 
speed over five minutes at a given height remains essentially constant for 
successive five minute periods; the maximum and minimum variations in 
the wind speed must occur within that five minute period. Assuming that 
the eddy pair is carried by the mean wind, the maximum horizontal length 
scale is just the mean wind times five minutes. 

The information known at this point allows a x and a z to be estimated. 

Next, the variational constant B of the eddy pair must be estimated. The 
value of B in equation (9) is varied until the fluctuation of the horizontal 
wind speed agrees with the maximum and minimum wind speeds recorded 
over a five minute period on the tower. With these three parameters 
estimated, the eddy pair model will provide the horizontal and vertical 
wind speed with range and height. 

Determination of Eddy Pair Parameters 

The meteorological data consisted of five minute averages with the 
maximum and minimum of the wind speed in that period. A direct 
calculation of the scale sizes of the eddy pairs can not be made since they 
typically passed the tower in less than five minutes. The procedure used to 
determine the eddy pair parameters outlined in the previous section is used 
for the experimental runs examined. 

The first experiment to be examined is Run 2.1 of January 11, 1985. 

The important constants are the mean wind speed, the horizontal and 
vertical wavenumbers, and the constant, B , for the eddy pair. The mean 
wind speed is calculated from the meteorological profiles of the 
experimental runs by performing a curve fit to equation (2). The procedure 
to determine the horizontal and vertical wavenumbers is to use equation (3) for 
calculating the height of the first inversion layer and using this height to calculate 
the vertical height of the eddy pair. The curve fit to equation (2) provides 
values for the roughness length and the friction velocity. The horizontal 
wind speed fluctuation, <r u , is determined from the hot wire measurements. 

Using the mean wind speed and incoming solar radiation, the 
Monin-Obukhov length can be estimated from table 2. 

For the experiment in question, the day was overcast with a light wind 
of 3.3 m/s. Using tables 1 and 2 for incoming solar radiation and a surface 
wind speed of 3.3 m/s, the Monin-Obukhov length, T mo , was estimated to 
be 20 meters. From analysis of the five minute wind speed measurements 


with height, the horizontal wind speed fluctuation was 0.40 m/s. For the 
experiments discussed, the roughness length was estimated from table 3 to 
be 0.05 meters. Using these parameters, the height of the first inversion 
layer was calculated to be 450 meters using equation (3). 

Using the condition that the eddy pair traverses past the tower within a 
five minute period, the maximum eddy pair size possible to traverse the 
field of propagation is 990 meters. If the dimensions of each eddy are 450 m, 
then the eddy pair has a horizontal length of 900 meters. This size is less 
than the maximum size constraint dictated by the five minute measurement 
period. Using equation (9), the horizontal and vertical wavenumbers ( a x 
and a z ) for the eddy pair are 0.125 m -1 and 0.025 m -1 . 

To determine the constant B in equation (9), the maximum and 
minimum wind speed fluctuations within a five minute segment with height 
are compared with the wind speed fluctuations predicted by the model. 

The parameter B is varied until the predicted wind speed variations fit 
those observed for a five minute segment. For the date in question, the 
value of B which best fit the data is 200 m 2 /s. 

The next experimental run was Run 1.1 of December 13, 1984. This day 
differed from January in that the mean wind speed and horizontal wind 
speed fluctuations were much greater. The mean wind speed was 6.3 m/s 
while the horizontal wind speed fluctuation was 1.0 m/s. Table 6.6 in 
Panofsky and Dutton 1 is used to determine the value of L mo . Using this 
table, the value of L mo is estimated to be on the order of 100 to 150 m, 
which gives a value for z, of 575 to 875 m. 

Results From the Eddy Pair Model 

Viewing the movement of the eddy pair on the scale of the geometry of 
the experiments, the variation of the sound speed in the atmosphere would 
appear to change slowly over the entire range of the experiment uniformly. 
Using a simple model of the wind speed in the atmosphere slowly varying 
from U\ to U 2 , the expected phase change can be calculated using 

A - « 2 ) (11) 

where R is the propagation distance, c 0 is the sound speed at temperature 
T, and / is the frequency of the signal. A comparison between the 
magnitude of the phase change for the simple model and the experimental 
results is shown in table 3. 

Conclusions 

Experimental acoustic phase data definitely displays two variational 
time scales. The short term time variations can be attributed to the 
presence of small scale turbulence present in the atmosphere. The small 
scale turbulence does not account for the longer time variations in phase. 


The large scale turbulence model is composed of pairs of vortices or 
eddies moving through the atmosphere at the mean wind speed. The scale 
parameters for the eddy pairs are determined from the available 
meteorological data composed of the maximum, minimum, and average 
wind speed over a five minute segment for four heights and meteorological 
theories of the behavior of the lower atmosphere. The constraint of the 
eddy pair moving through the field of propagation within five minutes is 
generally used as an upper bound for the dimensions of the eddy pair; 
however, it could be used as the size of the eddy pair if there is lack of 
available meteorological data. 

The results of the eddy pair model were used to examine the phase 
fluctuations of the acoustic wave using a simple phase model. The input 
parameters for the model were determined from analysis of the acoustical 
and meteorological data collected in the experiments. The magnitude of the 
phase variations predicted using this model was found to be in very good 
agreement with the experimental results. 
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Table 1 . Estimation of Turner Classes. 


Surface Wind 

Incoming Solar Radiation 

Speed (at 10m), m/s 

Strong Moderate 

Light 

<2 

1 1 

2 

CO 

i 

C\J 

1-2 2 

3 

3-5 

2 2-3 

3 

5-6 

3 3-4 

4 

6< 

3 4 

4 

lable 2. Estimation of L mo for Various Turner Classes. 


Turner Class 

" Lmo 


1 

8-12 m 


2 

12-20 m 


3 

20-60 m 


4 

>60 m 
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Table 3. Results from the Simple Phase Model for Run 2.1 of January. 
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Figure 1. Breakdown of The Lower Atmosphere. 
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Figure 4a. Time Fluctuations for High Wind Speed. 
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Figure 4b. Spectrum of Phase Fluctuations for High Wind Speed. 
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